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SECTION  I 


INTRODUCTION 

The  past  two  decades  have  witnessed  a  rapidly  expanding  interest  in 
characterization  of  high  speed  flows,  principally  in  aerodynamics.  This 
has  fostered  formulation  and  use  of  numerical  algorithm  constructions  for 
the  governing  partial  differential  equation  system.  Coupled  with  recent 
advances  in  computer  hardware/f irmware  capabilities,  it  appears  that  the 
age  of  realistic  computational  fluid  dynamics  (CFD)  analysis  may  be  nearer 
at  hand.  In  concert  with  evolution  of  the  three-color  laser  velocimeter, 
both  experimental  and  theoretical  fluid  dynamics  analysis  should  soon  enjoy 
a  greatly  expanded  resolution  capability. 

The  true  "workhorse"  of  the  CFD  aerodynamics  community  over  the  past 
decade  has  been  the  "MacCormack  algorithm"  and  variations  thereof.  For 
example,  fully  half  of  the  CFD  results  presented  at  the  recent  NASA  Lewis 
Research  Center  Workshop  [1]  were  generated  using  this  algorithm,  originally 
published  in  1968  [2].  This  is  due  principally  to  the  ultimate  simplicity 
of  this  explicit,  predictor-corrector  algorithm,  and  its  proven  track  record 
for  prediction  of  compressible,  supersonic , inviscid  flowfields.  In  the  split- 
operator  construction,  the  programming  requirements  are  elementary,  and  the 
resultant  code  runs  quite  economically.  The  basic  theoretical  formulation 
enjoys  continuing  refinement,  including  application  to  the  viscous  Navier- 
Stokes  equations,  cf.  [3,4]. 

The  numerical  simulation  of  viscous  or  turbulent  aerodynamic  flows 
places  considerable  additional  demands  on  a  CFD  algorithm.  The  added  dis¬ 
cretization  refinement,  required  to  resolve  wall  layers,  "stiffens"  the 
resulting  discretized  equation  system.  For  an  elementary  explicit  algorithm 
construction,  the  integration  absolute  stability  interval  is  bounded  by  a 
multiple  of  the  largest  eigenvalue  of  the  matrix,  yielding  the  requirement 
to  time  march  using  a  very  small  step  size.  For  this  reason  at  least,  and 
on  its  own  merit,  CFD  research  attention  has  turned  in  the  last  five  years 
to  development  of  implicit  Navier-Stokes  algorithms.  Simply  stated,  one 
trades  the  small  time-step  restriction  for  the  requirement  to  solve  matrix 
equation  systems.  As  a  result,  principal  attention  has  focused  on  matrix 
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factorization  procedures  that  reduce  largesparse  matrices  to  block-banded 
forms  to  permit  running  solutions  on  present  computer  systems. 

The  approximate-factorization,  implicit  finite  difference  (AFFD) 
Navier-Stokes  algorithms  of  Beam  and  Warming  [5],  and  Briley  and  McDonald 
[6],  exemplify  the  concept.  The  numerics  are  typically  second-order  accurate 
in  space,  first  order  in  (pseudo-)  time,  and  may  employ  both  "implicit"  and/ 
or  "explicit"  artificial  diffusion  to  control  perceived  instabilities.  In 
the  favored  non-iterative  "delta"  formulation,  the  matrix  solution  proce¬ 
dure  constitutes  acceptance  of  the  first  iterate  of  the  Newton  algorithm, 
using  an  approximation  to  the  true  matrix  system  Jacobian,  constructed 
upon  addition  of  the  identity  matrix. 

An  integral  ingredient  of  the  AFFD  algorithm  construction  is  definition 
and  use  of  the  "generalized  coordinates"  description.  Basically,  the  diver¬ 
gence  operator  in  the  Navier-Stokes  equation  set  is  transformed  into  scalar 
components  parallel  to  principal  coordinates  of  a  coordinate  transformation 
"regularizing"  the  boundary  3R  of  the  solution  domain  P.n  to  the  unit  square 
(or  cube,  n=3).  Thompson  and  coworkers  [7]  pioneered  the  transformation 
concept,  the  subject  of  which  has  become  rather  popular,  cf.  [8].  Steger 
and  Pulliam  [9]  first  reported  results  using  the  generalized  coordinates 
AFFD  algorithm  concept  for  a  two-dimensional  cascade  flow;  topical  results 
were  recently  summarized  by  Steger  [10]. 

In  this  same  time  period ,  serious  attention  has  become  focused  on  exami¬ 
nation  of  finite  element  concepts  applied  to  CFD.  In  its  classical  applica¬ 
tion,  the  "finite  element  method"  is  the  engineer's  utilization  of  the 
mathematician's  formulation  of  the  variational  boundary  value  problem. 
Specifically,  given  a  functional  describing  the  state  of  a  system,  extremize 
it  with  respect  to  all  eligible  states,  to  determine  the  constraint  set  of 
the  extremum.  The  numerical  solution  of  this  matrix  equation  system  yields 
the  engineering  parameters  of  interest.  Embedded  within  is  the  variational 
calculus,  and  the  formality  is  applicable  to  a  wide  range  of  problem  classes 
in  mechanics  and  physics,  cf.  [11,  12]. 

In  its  more  elementary  interpretation,  the  finite  element  (FE)  method 
returns  calculus  and  vector  field  theory  to  the  construction  of  discrete 
simulation  algorithms  for  any  branch  of  mechanics.  Of  necessity,  using  Tay- 
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lor  series  expansions,  one  must  always  be  able  to  verify  the  equivalent 
(finite  difference)  order-of-accuracy  for  elementary  derivatives  within 
the  governing  equation  system.  Specifically,  linear  (quadratic)  finite 
elements  yield,  respectively ,  second-  (  fourth-)  order  accurate  finite  dif¬ 
ference  representations  for  linear,  first  and  second  order  spatial  deriva¬ 
tives.  For  other  than  these  linear  spatial  derivatives,  the  finite  element 
construction  usually  produces  expressions  that  are  not  so  familiar.  However, 
upon  dissection,  they  can  always  be  related  to  an  appropriate  Taylor  series 
expansion  ,  hence  be  expressed  as  a  difference  recursion  relation.  The  im¬ 
portant  issue  is  that  the  theory  produces  the  discrete  analog  expressions, 
completely  independent  of  the  a  posteriori  ability  to  construct  an  equivalent 
difference  representation. 

This  time  frame  has  also  witnessed  emergence  of  the  "finite  volume" 

(FV)  algorithm,  cf.  Shang  [13].  A  FV  algorithm  is  constructed  using  direct 
integration  of  the  conservative-form,  governing  equation  system,  replacing 
derivatives  by  divided  (finite)  differences.  Fluxes  are  expressed  as  inte¬ 
grals  over  the  control  surfaces  of  a  FV,  and  scalar  fields  are  typically  cell- 
centered.  In  an  explicit  formulation,  the  FV  procedure  yields  the  geometric 
versatility  of  a  FE  discretization,  without  incurring  the  coordinate  trans- 
formation-induced  singularities  sometimes  associated  with  the  generalized 
coordinates  finite  difference  formulation.  An  implicit  FV  algorithm  would 
require  handling  of  large,  sparse  matrix  Jacobians,  unless  a  procedure 
were  invoked  to  yield  the  block-banded  structure. 

The  focus  of  this  research  project  has  been  to  address  some  of  the 
basic  issues  regarding  theoretical  and  practical  aspects  of  CFD  algorithm 
constructions,  in  particular  accuracy,  convergence  and  efficiency.  The  pre¬ 
vious  two  Interim  reports  [14,15]  document  certain  of  these  results,  and 
this  report  summarizes  the  principal  results  of  the  three-year  project. 

There  are  four  basic  requirements  that  must  be  met  in  the  construction  of 
any  algorithm  for  CFD,  see  Table  1.  Denoting  the  representative  partial 
differential  equation  of  the  Navier-Stokes  set  as  L(q(x,t)),  the  first  and 
most  obvious  requirement  is  that  an  approximation  must  be  specified.  The 
options  available  include  (at  least)  finite  element,  finite  difference  and/ 
or  finite  volume  concepts.  Since  the  approximation  cannot  be  the  exact 
solution,  a  statement  is  required  regarding  what  constraints  (if  any)  are 
placed  on  the  approximation  error. 


3 


Table  1 


Basic  Requirements  For  A  CFD  Algorithm  Construction 


Category 

Options 

1.  Approximation 

Finite  Element  (FE) 

Finite  Difference  (FD) 

Finite  Volume  (FV) 

2.  Approximation  Error 

Orthogonality  (FE) 

Taylor  Series  (FD,  FV) 

3.  Dispersion  Error 

Multi -pole  Expansion  (FE) 

Artificial  Diffusion  (FD,  FV) 

4.  Matrix  Solution  Jacobian 

Tensor  Matrix  Products  (FE) 

Approximate  Factorization  (FD) 

The  FD  and  FV  constructions  differ  from  FE  concepts  in  this  spe¬ 
cific  regard,  in  that  categories  1  and  2  are  grouped  together  in  replacing 
L(q (x  ,t) )  with  difference  quotients  yielding  L (Q (Ax  . ,At) )  =  0.  In  the  FE 

J 

construction  of  category  1,  the  approximation  is  made  to  the  solution  rather 
than  L(q(x,t)).  A  vast  reservoir  exists,  from  which  one  can  select  particular 
members  that  endow  the  solution  approximation  q^(x,t)  with  specific  proper¬ 
ties,  eg.,  completeness.  Then,  L(q^(x,t))  f  0  is  the  error  in  the  governing 
differential  equation.  The  category  2  specification  simply  requires  this 
error  to  be  orthogonal  (perpendicular)  to  the  function  space  from  which  q*1 
was  formed.  Calculus  and  vector  field  operations  permit  an  exact  evaluation 
of  this  constraint,  yielding  the  matrix  equivalent  statement  L({Q(t)})  =  {0}. 

For  a  viscous-dominated  (small  Reynolds  number)  flow,  the  category  3 
issue  of  dispersion  error  need  not  be  addressed.  However,  most  aerodynamic 
flows  are  characterized  by  a  large  Reynolds  number,  yielding  the  governing- 
Navier-Stokes  equations  elliptic  only  in  small  regions  of  Rn  directly  adja¬ 
cent  to  aerodynamic  surfaces.  Elsewhere,  the  equation  set  is  dominantly 
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hyperbolic  and  non-linear.  The  principal  error  mechanism  operating  in 
discrete  approximations  tbdiyperbolicequations lijsdispersion  error,  the  ten¬ 
dency  of  the  Fourier  components  constituting  information  packets  to  possess 
distinct  phase  velocities  differing  from  the  group  velocity.  The  Navier- 
Stokes  non-linearity  aggravates  this  character,  yielding  a  divergence  that 
can  eventually  "blow-up"  the  solution.  The  algorithm  requirement  is  to 
control  (diffuse)  this  error  mechanism.  The  standard  FD(FV)  practice  is 
to  explicitly  add  a  diffusion  term,  multiplied  by  a  coefficient  of  "arti¬ 
ficial  viscosity,"  to  the  governing  differential  equation  L(q(x,t)).  The 
concept,  first  formalized  by  von  Neumann  and  Richtmyer  [16],  is  broadly 
used.  The  FE  basis  for  this  error  control  is  recognition  that  the  approxi- 
mation  error  L(qr  (x,t)  will  be  highly  non-smooth  in  regions  where  dispersion 
error  is  excessive.  Interpreting  L(q^(x,t))  as  a  (hyper-)  surface,  regions 
of  large  error  are  characterized  by  large  changes  in  surface  amplitude  and 
direction  of  the  normal.  Noting-that  VL(qh(x,t))  characterizes  both  the 
magnitude  and  direction  of  the  non-smoothness,  the  FE  concept  requires  this 
error  also  be  orthogonal  to  the  approximation  sub-space,  subject  to  a  con¬ 
straint  set|.  This  of  course  yields  an  "artificial  viscosity"  term,  as 
well  as  some  additional  terms,  and  it  is  interesting  to  compare  performance 
of  the  FD  and  FE  formalisms. 

The  completion  of  categories  1-3  yields  a  matrix  equation  system 
{ F} ,  that  is  non-linear  unless  explicit  time-integration  has  been  specified. 
In  all  other  instances,  including  a  direct  steady-state  algorithm,  and  as  a 
function  of  the  dimension  n  of  Rn  and  the  number  of  dependent  variables.de- 
fined  in  L(q),  the  Jacobian  of  { F}  is  a  large  albeit  sparse  matrix.  The 
direct  matrix  solution  using  this  Jacobian  is  usually  ill-advised,  based 
strictly  on  the  operation  count,  core  storage  and  CPU  requirements.  One  FD 
approach  is  to  construct  an  approximate  factorization,  as  discussed  earlier. 
One  FE  approach,  presented  and  discussed  herein,  is  to  construct  the  tensor 
(outer)  matrix  product  of  the  Jacobian.  Of  course,  since  category  4  is 
strictly  a  problem  in  linear  algebra,  there  exists  a  wide  range  of  alterna¬ 
tive  procedures  including  SOR,  SLOR,  and  ADI  to  name  a  few.  The  selection 
of  any  of  these  can  be  a  (the)  prime  determining  factor  in  the  cost  of  ob¬ 
taining  the  solution.  The  solution  cost  is  typically  much  less  affected  by 
choices  made  in  categories  1-3. 
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Before  proceeding  into  the  detailed  technical  analysis,  an  amplifi¬ 
cation  on  category  3  is  appropriate.  The  concept  of  "artificial  viscosity" 
has  been  rediscovered  innumerable  times,  under  such  euphemisms  as  upwinding, 
upstream  differencing ,  donor  cel  1 ,  streamline  upwinding,  Petrov-Galerkin,  non¬ 
equal  weighting  (and  multi-pole  expansion).  In  every  instance,  terms  are 
added  to  smooth,  the  solution.  There  has  also  been  a  tendency  to  prescribe 
a  dose  of  artificial  viscosity  to  smooth  error  introduced  by  inadequate  grid 
resolution,  distorted  grids,  and  inconsistent  formulations,  or  to  "remedy" 
the  overall  poor  phase  accuracy  of  the  specific  algorithm.  Strictly 
speaking,  the  use  of  fictitious  diffusion  to  control  other  than  non-linear 
induced  instability  constitutes  an  additional  error  source  not  covered  in 
Table  1.  The  "suppression  of  /the  wiggles",  cf.  [17],  really  constitutes 
elimination  of  information  about  the  approximate  solution,  specifically  its 
inadequacy.  This  facet  of  CFD  certainly  deserves  close  attention,  in  the 
context  of  a  missing  category  in  Table  1. 


6 


SECTION  II 


PROBLEM  STATEMENT 

The  partial  differential  equation  set  governing  three-dimensional 
aerodynamic  flows  is  the  familiar  and  very  non-linear  Navier-Stokes  system. 
In  non-dimensional  conservation  form,  using  Cartesian  tensor  summation 
notation,  the  equation  set  governing  flow  of  a  compressible,  viscous,  heat- 
conducting  fluid  is 


=  IMrCY3  = 0 

«J 


L(pui ) 


3(pUj)  +3 


3t 


3x. 


u.  pUi  +  P6.  .  - 


=  0 


L(pe)  * 


+  ujp  '  aijui  “  qj]  =0 


(1) 

(2) 

(3) 


In  equations  1-3,  p  is  density,  pu^  =  nr  is  the  momentum  vector,  p  is  pres¬ 
sure,  and  e  is  mass  specific  total  energy.  The  Stokes  viscous  stress  tensor 
and  heat  flux (Vector  q j ,  in  terms  of  specific' internal  energy  e,  are 


o 


_  u 

3U^ 

ij  "Re 

3Xd 

v  j^k 

3Re  3xk 


(4) 


-K 


(5) 


e  =  e 


(6) 


Assuming  a  polytropic  gas,  p  =  (y-l)pe,  and  the  equation  of  state  is 


L(p)  =  p  -  (Y-i)[Pe  -  i-pUjUj]  =  0  (7) 

Finally,  y  is  the  absolute  viscosity,  k  is  the  coefficient  of  heat  conduc¬ 
tivity,  6 • •  is  the  Kronecker  delta,  and  Re  is  the  reference  Reynolds  number. 

•  J 
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The  Euler  equations  are  contained  within  this  description  by  the 
specification  that  equations  4  and  5  vanish  identically.  The  form  of  equa¬ 
tions  1-3  is  also  representati ve  of  the  mass-weighted  formulation  of  the  time- 
averaged  Navier-Stokes  equations  for  a  turbulent  flow  [18].  The  dependent 
variables  are  interpreted  as  descriptors  of  the  mean  flow,  and  ckj  is 
generalized  to  include  nonvanishing  correlations  of  sub-grid  scale  phenomena. 
In  this  instance,  the  total  stress  tensor  becomes  of  the  form. 


pUiUj 


(8) 


where  -puiu^  is  the  dynamic  Reynolds  stress  tensor,  and  a*,  denotes  the 
time  averaged  form  for  equation  4. 

Procedures  for  closing  equation  8  include  solution  of  the  Reynolds 
stress  transport  equations,  cf.  [18,19],  tensor  field  constitutive  theory 
[20],  as  well  as  elementary  mixing  length  theory  [18].  For  the  first 
closure,  a--  will  not  be  an  explicit  function  of  m . ;  for  the  latter  two,  it 

I  J  * 

becomes  a  strongly  non-linear  function  of  m^ .  The  field  theory  approach 
encompasses  the  latter,  in  the  form 


u-u: 
i  J 


ci“u 


riki 
v  e 


r2Ji- 
v  e 


au}. 
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3X 


(9) 


The  various  correlation  coefficients  are  defined  [20],  and  the  turbulence 
kinetic  energy  and  isotropic  dissipation  function  are,  respectively 


k  = 


e  = 


1  u^ur 

2  ii 


2v 

9u  1 

3 

b 

\jk 


(10) 

(11) 


The  variables  k  and  e  are  solutions  to  the  corresponding  transport  differen¬ 
tial  equations  [19], 
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(12) 


(13) 


The  correlation  coefficients  and  have  been  determined  from  analysis 
and  experiment  [19]. 

The  boundary  conditions  for  the  equation  set  1-13  are  a  general  mix¬ 
ture  of  Dirichlet  and  Neumann  specifications.  On  an  inflow  boundary  segment, 
p,  pu-,  pe,  k  and  e  are  typically  specified  and  p  is  determined  from  equa¬ 
tion  6.  For  an  inviscid  streamline  (Euler  equations  on  a  farfield  boundary), 
the  normal  derivative  of  the  scalar  variables  vanish,  and  pu..is  con¬ 
strained  to  be  parallel.  For  a  no-slip  aerodynamic  surface,  pu^ ,  k  and  e 
vanish  identically,  and  the  normal  derivative  of  0  is  constrained  (adiabatic 
or  cooled  wall).  At  an  outflow  boundary,  typically  p  is  specified  and  the 
normal  derivative  of  the  remaining  dependent  variables  is  constrained. 
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SECTION  III 


FINITE  ELEMENT  SOLUTION  ALGORITHM 

1 .  Theoretical  Statement 

Equations  1-13  constitute  a  fairly  general  definition  of  the  three- 
dimensional  Navier-Stokes  equations.  In  the  limit  Re  =>  equation  4,  they 
reduce  to  the  Euler  equations  governing  an  inviscid  rotational  flow.  Spec¬ 
ifying  m.  s  pu.j,  g  =  pe,  define  the  vector-valued  dependent  variable  set  as 


qa(xi,t)  5  {qa>  =  {p,mi >g,p,a. j ,qj ,  pk,pe) 


(14) 


Equations  1-5,  12,  13  constitute  an  initial-value,  elliptic  boundary  value 
description  of  the  general  form 


L(qa) 


+  s  =  0 

a 


(15) 


In  equation  15,  f^q^)  and  sa(q^)  are  specific  non-linear  functions  of 
their  argument,  as  determined  by  inspection.  The  remaining  algebraic  and 
partial  differential  equations  7-9  are  of  the  form 


L(qa>  =  qa  +  fa(qB)  =  0  (16) 

where  the  f  are  also  determined  by  inspection.  Note  that  equation  16  is 
simply  a  special  form  of  equation  15. 

The  n-dimensional  partial  differential  equation  system  15-16  is 
defined  on  the  Euclidean  space  Rn,  spanned  by  the  x  coordinate  system  with 
scalar  components  x. ,  1  <  i  <  n.  The  solution  domain  ft  is  defined  as  the 
product  of  Rn  and  t,  for  all  elements  of  x  belonging  to  Rn  and  all  elements 
of  t  belonging  to  the  open  interval  measured  from  t  ,  i.e.. 


The  boundary  3fl  of  the  solution  domain  is  the  product  of  the  boundary 
3R  of  Rn ,  spanned  by  t,  and  t ,  i .e. ,  3^  =  3R  x  t.  Thereupon,  the  generalized 
form  for  the  differential  boundary  constraint  is 


H(%)  =  a?q„  + 


a  3 


'a 


a 


2  3x. 


q  n. 

j 


+  <  *  ° 


(18) 


In  equation  18,  the  a^  are  specified  coefficients  and  n.  is  the  outwards 

pointing  unit  normal  vector.  Finally,  an  initial  distribution  for  q  on 

n  a 

=  R  xtQ  is  required. 


■  •£<*>  <19> 

The  functional  requirement  of  the  (any)  numerical  solution  algorithm 
for  equations  14-19  is  to  construct  a  suitable  "approximation"  to  the 
dependent  variable  set  q  (x,t),  and  to  render  the  error  in  this  approximation 
extremum  in  some  norm,  recall  Table  1.  A  finite  element  algorithm  construc¬ 
tion  is  simply  the  formalized  application  of  this  basic  requirement,  as 
discussed.  The  semi-discrete  approximation  q  (x,t)  to  the  unknown  exact 

solution  q  (x,t)  is  constructed  from  members  of  a  convenient  f inite-dimen- 

a 

sional  subspace  of  H0(R),  the  Hilbert  space  of  all  functions  possessing 
square  integrable  first  derivatives  and  satisfying  the  boundary  conditions, 
equation  18.  While  extremely  flexible  in  theory,  the  typical  practice  is 
to  employ  elementary  polynomials,  truncated  at  degree  k,  and  defined  on 
disjoint  interior  subdomains  R^,  the  union  of  which  forms  the  discretization 
of  Rn  =  UR".  Hence, 

.  M 

q  (x,t)  «  qjj(x,t)  =  E  q5(x,t)  (20) 

a  a  e=l 

and 


q®(x,t)}  5  (Nk(x)}T{QI(t)}e  (21) 

In  equation  21,  the  semi -discrete  free  index  "I"  denotes  q^1  evaluated  at 

n  a 

the  nodal  coordinates  of  the  discretization  URg  at  anytime  t.  The  sub-  or 
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superscript  e  indicates  pertaining  to  the  e—  finite  element  domain, 

=  Rext*  The  alements  the  row  matrix  (Nk(x)}T  are  polynomials  written 
on  Xj,  1  <  j  <  n,  complete  to  degree  k  and  constructed  to  form  a  cardinal 
basis  [21]. 

The  fundamental  requirement  of  any  numerical  algorithm  is  to  render 

the  (semi -discrete)  approximation  error  minimum  in  some  norm.  For  smooth 

solutions,  this  is  accomplished  within  a  finite  element  context  by  requiring 

h  h 

the  error  in  equations  15,  16  and  18,  ie.,  L(qn)  and  £(qn),  to  be  orthogonal 

l.  ot  ot 

to  the  function  space  used  to  construct  q  .  This  basic  requirement  may  be 
augmented,  in  the  sense  of  a  multi -pole  expansion,  to  enforce  further  con¬ 
straints  on  the  algorithm  construction.  For  example,  for  high  Mach  number 
flows  wherein  non-smooth  solutions  may  exist,  the  dissipation  mechanism  for 
control  of  non-linearly  induced  instabilities  is  generated  by  also  requiring 
the  error  (vector)  gradient  VL(q^)  to  be  orthogonal  to  {N^} .  As  a  second 
example,  for  certain  steady,  subsonic  flows,  the  solution  to  the  momentum 
equation  2  can  be  constrained  [20]  by  the  requirement  that  the  velocity  field 
be  divergence  free,  equation  1.  Each  of  these  constitutes  a  term  in  a  multi - 
pole  expansion  of  the  basic  error  orthogonal ization.  Defining  the  (Lagrange) 
multiplier  set  ]§.  permits  adjoining  these  linearly  independent  constraints, 
yielding  the  finite  element  numerical  solution  algorithm  statement  for  the 
compressible  Navier-Stokes  equations  as 


{Nk}L(q>^  +  fc 


9R 


{NkH(q5)dx  +|2- 


'a 


{Nk}VL(q>;  =-  {0} 


(22) 


Upon  definition  of  k  in  equation  21,  equation  22  is  a  mixed  system  of 
ordinary  differential  and  algebraic  equations,  written  on  t,  of  the  form 


[C]{QIK  +  [U]{QI>  +  [FI  J  ]  {QJ }  +  {SI}  =  {0}  (23) 


A  one-to-one  correspondence  of  terms  in  equations  23  and  15-16  is  inferred, 
with  the  matrices  in  equation  23  augmented  to  contain  the  various  additional 
terms  introduced  through  6^  t  0  in  equation  22.  (Detailed  expansions  are 
presented  in  a  later  section.)  Since  the  [C]  matrix  is  non-diagonal,  the 
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form  of  equation  23  suggests  use  of  an  implicit  integration  algorithm.  The 
efficient  0-implicit  finite  difference  algorithm,  where  0  =  i  yields  the 
trapezoidal  rule,  is 


{FI}  =  {QI} 


j+1 


-  «1}j  - 


At 


e{Qi}j+1+  (1-0)  {qi  }  T 


=  0  (24) 


Equation  23  provides  the  definition  of  the  derivatives  {QIK  Upon  substi¬ 
tution,  and  proceeding  through  the  algebra  [14],  equation  24  becomes  a  non¬ 
linear  algebraic  equation  system  written  on  {QI ) j+x *  The  Newton  iteration 
algorithm  solution  statement  is 


[J  (FI )  ]  j+1{6QI  >  j+i  =  -{FI}[?+i  (25) 

The  dependent  variable  in  equation  25  is  the  iteration  vector  {6QI } ,  related 
to  the  solution  in  the  conventional  manner. 


«!>$  5  {QI}?+1  +  MM 


P+1 

j+1 


(26) 


The  Jacobian  is  defined  as, 


=trtOT 


(27) 


The  detailed  construction  is  direct,  upon  expansion  of  equations  22-24 
using  a  hypermatrix  formulation,  to  be  discussed. 

2 .  Generalized  Coordinates 

A  principal  requirement  in  computational  aerodynamics  is  to  accurately 
interpolate  domain  boundary  geometries  that  are  non-regular  and  perhaps 
non-smooth.  The  term  "generalized  coordinates"  has  gained  acceptance  to 
describe  the  algorithm  statement  appropriate  for  use  with  a  regularizing, 
boundary-fitted  coordinate  transformation.  Many  procedures  are  available 
to  construct  such  transformations  [8],  including  numerical  solution  of 
elliptic  and  hyperbolic  partial  differential  equations  and  algebraic  methods. 
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The  output  of  any  of  these  procedures  constitutes  definition  of  coor¬ 
dinate  triples  (pairs)  on  R3  (R2),  that  define  the  intersections  of  3R  on 
Q  called  nodes.  The  computational  requirement  is  to  construct  a  local 

coordinate  transformation,  using  these  data,  that  maps  x^R^  to  a  regularized 
domain  spanned  by  an  orthonormal  coordinate  system  rj..  Figure  1  illustrates 

J 

the  concept,  where  (•)  and  (x)  depict  the  nodal  coordinates  in  both  spaces. 
The  required  coordinate  transformation  is. 


x 


i 


xi(V 


(28) 


The  explicit  form  for  equation  28  is  constructed  using  elementary  interpola¬ 
tion  concepts  and  the  cardinal  basis  >.  Since  equation  21  defines  any 
approximation,  the  form  for  equation  28  is 


x.j  =  {Nk(q))T{XI}e,  x1-eRg 


(29) 


a)  Two-Dimensional 
Domain 


b)  Three-Dimensional 
Domain 


Figure  1:  Biquadratic  Cardinal  Basis  Coordinate  Transformation 
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e> 


The  elements  of  the  column  matrix  {XI >e  are  the  coordinates  of  the  nodes  of 
Rg,  1  <  (i,I)  <  n,  hence  UR^,  the  global  discretization. 

Assuming  the  discretization  of  Rn  adequately  represented  by  plane¬ 
faced  hexahedra  (or  straight-sided  quadrilaterals,  for  R2) ,  the  k=l  basis 
{Ni(o)>  yields  the  exact  representation.  For  example,  on  R2, 


{Ni(n)>  =  \ 


jl-m)  (1-02)] 
J  (1+m)  (l-m) [ 

I  (1+m)  (1+02) f 
(l-m)(i+m) 


(30) 


A  discretization  on  R2  using  curved-sided  quadrilaterals  must  contain  at 
least  eight  nodes  per  element,  and 


{Na(n)> 


1 

4 


(l-m)  (l-m) (-01-O2-I) 
(1+m )  (l-f|2 )  (  ni-r|2-l) 
(1+m) (1+02) (  m+m-l) 
J  (l-m)  (1+m)  (-m+m-l) 
]  (l-m)  (l-m) 
(l+m)(l-m) 

(1+m)  (1+02) 

(l-m)  (1+m) 


(31) 


For  discretizations  of  R3  employing  planar-faced  hexahedra,  the  corresponding 
cardinal  basis  is 


(Ni(r[)}  =  ^ 


(1-Pl)  (I-O2)  ( 1  -T|  3  ) 

(1+m)  (l-m)  (l-m) 
(1+m) (1+T12 )  (l-m) 
J  ( 1-rii )  (1+m)  ( l-r)3 ) 

1  (1-Tll)  (1-T12)  (1+03) 
(1+m)  ( 1  -r)2 )  (1+m ) 
(1+m)  (1+T12 )  (1+m ) 
(l-m)  (I+02 )  (1+03) 


> 


j 


(32) 


Returning  to  the  issue,  the  generalized-coordinates  algorithm  requirement 
is  transformation  of  the  divergence  operator  in  equation  15,  i.e., 


9 


9x, 


!Hi3_ 

3xi  3rij 
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(33) 


The  elements  of  the  inverse  Jacobian  J_1  =  [3n./3x.]  are, 

*  J 

=  J"1  -  [transformed  cofactor  of  J]  (34) 

J  is  the  Jacobian  of  the  forward  transformation,  equation  28, 


[J]  = 

which  is  directly  evaluated  in  terms  of  the  nodal  coordinates  of  using 
equation  29.  The  differential  element  for  equation  22  is 

dx  =  det[J]dri  (36) 


9xi 


9T\j 


(35) 


As  an  example,  consider  equation  22  written  for  the  momentum  equation  2. 
Using  a  Green-Gauss  form  of  the  divergence  theorem,  the  first  constraint 
term  of  equation  22  becomes 


h  3pu,M 

n(N}L(puJ)dx  =  n{N}  — det[J]dn 
R  ^  R 


+  0  {N} 

3R 


uV  +phS..  -  0h 

J\T  1  H1J 


igj  njdet[J]d1 


3{N} 

K] 

L8xjJ 

[ujpud  +  pd6i j  -  o'jj  det[J]drl  (37) 


Define  the  contravariant  components  of  the  convection  velocity  as. 


uj!  5  det[J] 


8Xj. 


“S 


(38) 


with  scalar  components  parallel  to  the  coordinate  system.  Using  the  alge¬ 
braic  transformation  equation  29,  det[J]  cancels  yielding 
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ujj  =  [Cof.  J]  uh. 


(39) 


where  [Cof.  J]  is  the  transformed  co-factor  matrix  of  [J].  Using  equation 
21  to  interpolate  the  nodal  distribution  of  ujj  on  r",  equation  37  becomes 


n{N}L{puV)dx  =  S€ 


n{DET}^{NHN}{N}T{RHOUI}gdn 


{UBARK}J(N}  {NHN}T{RHOUILd?i 


3n 


e 

RnLl-m'UJel"J  3nk 


{ETAKI }T{N}  (N}(N}T{P») 


„ { ET AKJ >X { N }  4—  {N}(N}T{SIGIJ}  ctf, 


3n 


{UBARK}T{N}(N}{N}T(RHOUI>  + 

3FL  0  3R 


{DET}^{N}{N}{N}T({P}e6IK  -  (SIGIK>e) 


rr- 

.nkdn 

(40) 


In  equation  40,  since  the  (N.)  for  qh  are  defined  locally  on  r4  the 
limits  for  the  defined  integrals  are  ±1.  S0  is  the  matrix  operator  projecting 
computed  element  contributions  to  the  corresponding  global  matrices,  equa¬ 
tion  22.  The  determinant  of  J  is  interpolated  on  the  element  R^  using  (N^} 
and  the  nodal  values.  Similarly,  Sn^/Sx^  is  interpolated  as  (ETAKI . 

The  e-subscripted  terms  in  equation  40  are  independent  ofr)k,  hence 
can  bei  extracted  from  the  integrand  leaving  only  products  of  the  polynomials 
and  derivatives  of  (N^Kcf.  equations  30-32.  These  are  directly  evaluated 
using  numerical  quadrature.  Hence,  the  finite  element  algorithm  statement 
for  the  first  term  in  equation  22,  as  formed  for  equation  2,  is 


’  (N}L(pu?)dx  =  Se 


I  (DET }g[M3000 ] (RHOUI }'  -  (UBARK}^[M30K0](RH0UI }p 

-  (ETAKI }g[M30K0]{P}e  +  (ETAKL}^[M30K0](SIGIL}e 
+  nJ(UBARK}*[B3000](RH0UI}  +  (DET}g[B3000](P}e6IK 


-  (DET}g[B3000](SIGIK}e 


(41) 
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In  equation  41,  the  indices  K  and  L  obey  the  tensor  summation  rule,  I 

L. 

is  the  free  index  for  pu?,  Sg  is  the  assembly  operator,  and  [M30K0]  is  the 
element-independent  hypermatrix  equivalent  of  a/Sri^  integrated  over  R^.  The 
matrix  [83000]  results  from  interpolation  integrals  on  the  element-boundary 
intersection  8Rg  A  3R  with  outward  normal  n^.  For  the  coordinate  transforma¬ 
tion,  (DET}g  is  the  nodal  distribution  of  det[J]  and  {ETAKI}g  and  (ETAKL}g 
are  corresponding  nodal  distributions  of  components  of  J-1  on  Rn.  Within 
the  generalized  coordinate  framework,  therefore,  the  grid  and  metric  data 
required  for  a  numerical  simulation  are  the  nodal  distributions  of 
J-1  =  [Sru/Sx-],  and  det[J]  =  det[3x-/3rw]>  These  data  are  strictly  a 

K  J  IK 

function  of  the  nodal  coordinate  distribution  of  the  discretization  URg, 
see  equation  28-29,  and  may  be  provided  using  any  coordinate  transformation 
procedure,  cf.  [8]. 

3.  Tensor  Matrix  Product  Jacobian 

For  efficiency,  recall  category  4  of  Table  1,  it  isdesirable  tocon- 
struct  a  suitable  approximation  to  the  equation  system  Jacobian  defined  in 
equation  25,  see  equation  27.  The  matrix  tensor  product  construction  is 
facilitated  by  using  the  tensor  product  cardinal  basis  function  set  (N^n)} 
spanning  quadrilateral  and  hexahedra  element  domains  Rg  and  R^  respectively. 
The  Jacobian  matrix  [J ( FI ) ] ,  equation  27,  is  then  replaced  by  the  tensor 
(outer)  product  construction  [22]. 

[J ( FI ) ]  =>  [Jx]  ®  [J2]  ®  [J3]  (42) 

Each  component  [J^]  is  constructed  from  its  definition,  equation  27,  assuming 
interpolation  and  differentiation  are  one-dimensional.  Using  equation  42, 
the  matrix  statement  equation  25  becomes 

[Jxl  ®  [02]  ®  [03]  {6QI>5+}  =  -  (FI>5+1  (43) 


Define 


[Jz]  <8>  [03]{6QI>5+J  e  {P1}?+1 
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[J,HSQI>$  =  {p2}^J 


(44) 


Then,  the  operation  defined  by  equation  43  is  the  sequence. 


[JiHPUjjt}  = 

-  {F»j+l 

CJ£]{P2}^j  = 

C'Js]{6QI}j+J  = 

<P2>$ 

(45) 

Obviously,  other  permutations  of  the  index  structure  for  [J^]  could 
be  utilized.  The  key  aspect  is  replacement  of  the  very  large  (albeit  sparse) 
Jacobian  matrix  [J],  with  a-block-structured  matrices  [J^].  The  principal 
attributes  are  several  orders  of  magnitude  reduction  in  central  memory 
requirements  for  the  Jacobian,  and  significantly  reduced  CPU  for  the  LU 
decomposition  and  back  substitution  solution  steps.  The  principal  detrac¬ 
tion  is  potential  degradation  of  the  quadratic  convergence  rate  for  the 
Newton  iteration,  equation  25.  This  procedure  in  no  way  affects  the  forma¬ 
tion  of  {FI},  equation  24,  wherein  lies  the  accuracy  features  intrinsic  to 
the  finite  element  algorithm  statement.  Compromises  in  the  evaluation  of 
{FI}  will  invariably  produce  inferior  results  for  the  Navier-Stokes  equations. 

The  construction  of  the  [J^]  is  straightforward.  For  example, 
for  the  initial -value  term  in  equation  41,  the  corresponding  term  expressing 
self-coupling  in  the  Jacobian  is 


=  [JQQD  =  Se  {DET }g [M3000]6 


(46) 


where  6j j  is  the  discrete  index  Kronecker  delta.  Referring  to  equation  40 
for  the  definition,  the  elemental  operation  in  forming  equation  46  is 


[JQQle  =  {DET}J[M3000]  = 


n{DET}^{Nk(^)}{Nk(n)}{Nk(n)}Tdn 

R_ 


(47) 
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Assume  for  simplicity  the  most  elementary  case,  i.e.,  k=l,  n=2  and  x  =  r\, 
i.e.,  the  identity  transformation.  Defining  M  =  B  for  n=2,  equation  47 
becomes 


{DET}g[B3000]  =  Ae[B200]  5 


JR 


(Ni(x)HNi(x)}Tdx 

2 

e 


(48) 


Assuming  the  rectangular  element  domain  described  by  measures  £  and  w, 
the  evaluation  of  equation  48  yields 


[JQQ3e  =  Ae[B200]  =  £w[B200]  =  || 


4  2  12 
4  2  1 
4  2 
Jsym)  4 


(49) 


The  tensor  product  construction  for  this  matrix  involves  evaluation 
of  equation  48  constrained  to  one-dimension.  Denoting  M  =  A  for  n=l. 


and 


[JQQ  ]  =  A  [A200]  =  {N1(x1)}{Ni(Xi)>  dxx 

n  e  n  JRi 


(50) 


[jqqne  =  f 


2  1 
1  2 


[JQQ2]e  =  f 


2  1 
1  2 


(51) 


assuming  Ai  =  £  and  A2  =  w.  Accounting  for  entry  locations  in  [JQQ], 
using  the  assembly  operator  Se,  cf.  [23],  it  is  readily  verified  that 

[JQQ]  =  Se 


[JQQi]e  ^ 


[jqq2]£ 


(52) 


By  direct  extension,  then,  the  tensor  matrix  product  form  of  equa¬ 
tion  46  on  R3  is 


(DET }g[M3000] 


=>  S  „  [{  D  ET 1 }"[ [A3000 ]  ® 


e[_ - eu 

[DET2}g[A3000]  ®  {DET3}T[A3000]  1 


(53) 
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In  equation  53,  the  index  1  <  K  <  3  on  {DETK}  denotes  the  scalar  component 
of  the  multi-dimensional  measure  parallel  to  r^. 

The  second  term  in  the  typical  Jacobian,  as  developed  in  Section  III. 4, 
is  the  form. 


,  v^{ETAKJ}g[M40K00]{DET}e 


Here  results  from  definition  of  the  specific  form  for  p2,  equation  22, 
and  a  denotes  the  specific  dependent  variable.  The  indices  1  <  (K,J)  <  n 
are  discrete  tensor  summation  indices  with  range  1  s  (K,  J)  <  n. 

[M40K00]  is  a  hypermatrix  of  degree  two,  ie.,  it  possesses  elements  which 
are  themselves  square  matrices,  which  represents  the  integration  over  of 
four  cardinal  basis  {^(n)},  one  of  which  is  differentiated  with  scalar  com¬ 
ponents  parallel  to  r^.  The  tensor  matrix  product  equivalent  is 


S 


e 


{ETAKJ  [M40K00]  {DET>e 


=>  S, 


{ ETA1 J }g [A40 100 ] {DET 1 }0  ®  <ETA2J}^ 


[A40l00]{DET2}e  ® 


{ETA3J}g[A40100]{DET}e 


(54) 


In  equation  54,  J  remains  a  discrete  summation  index  with  range  1  <  J  <  n. 

4.  Algorithm  Matrix  Statement 

As  developed  in  the  next  Section,  "§2  =  {DET)evJj  is  the  functional 
form  for  the  multi-pole  term  in  equation  22.  With  these  identifications, 
the  specific  matrix  forms  of  the  algorithm  statement  {FI}  =  {0},  equation 
24,  and  the  tensor  product  Jacobian  [J^],  equations  27  and  42,  can  be  ex¬ 
pressed.  In  all  cases,  the  global  matrix  statement  is  obtained  using  SQ 
on  the  sequence  of  elemental  definitions,  e.g., 

{FI}  ;  Se QfI)J  (55) 
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The  elemental  expression  { FI >g  is  always  obtained  using  the  element-inde¬ 
pendent  master  hypermatrices  [M. ..],  contracted  with  element-dependent 
matrices.  In  the  following  expressions,  subscripts  "e"  are  omitted  in 
the  expansions  for  clarity. 

Denote  the  members  of  the  discrete  dependent  variable  set  as  {QI}^  = 
{R,  M(I),  G,  P,  S (I ,J ) ,  Q(J),  RK,  RE},  see  equation  14.  The  matrix  algo¬ 
rithm  statements  { FI }g ,  equation  55,  for  arbitrary  basis  {N^}  and  dimension 
n,  are: 


(FR}e  = 


{DET}T[M3000]  +  vrj { ET AKJ }T [ M40K00 ] {D ET }  {R}j+1 


At 

~2~ 


-{ETAKI}T[M30K0]{MI}+  ^J{ETAKJ},[M40KL0]{UBARL}.{R}|j+1  .  (56) 


{ FMI }, 


(DET)T  [M3000]  +  Vjj{ETAKJ}T[M40K00]{DET} 


+  ^-MUBARK}T[M30K0]{MI}  -  {ETAKI}T[M30K0]{P} 


+  {ETAKJ}T[M30K0]{SIGIJ} 


Vjj  {ETAKJ  } 1  [M40KL0]  (UBARLHMI } 


j+U 


(57) 


{FG}e  = 


_,{gk 

-i  j+i 


{D  ET }T[M3000]  +  vqj(FTAKJ}T[M40K00]{DET} 

+  ^ -{UBARK}T[M30K0]({'G}+{P})  +  {ETAKJ}T[M30K0]{QJ} 
+  {ETAKJ}T[M40K00]{UL}{SIGLJ} 


+  V, 


GJ 


{ETAKJ}T[M40KL0]{UBARL}{G7|j+1  . 


(58) 
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{ FP>e  =  {DET }T[M3'G00]{P}  -  (y-1)  {DET}T[M3000]{6} 

+  HDET}T[M40000]{UJHMJ}  .+1  (59) 

(FSIJ)e  =  {DET}T[M3000]{SIGIJ}--^  (ETAKJ }T[M3OK0] {UI } 

f  +  {ETAKI}T[M30K0](UJ}  -  |(Sij{ETAKL}T[M30K0]{UL}Jj.+1  (60) 

{FQJ)e  =  (DET }T[M3000] (QJ }  -  k  {ETAKI{T[M30K0]{GSR}  -  {UJUJ}  ,+1  (61) 

{FRK)e  =  (DET}T[M3000]  +  v*j {ETAKJ }T[M40K00] {DET}  {RK}J+1 

+  ^p4UBARK}T[M30KO]{RK}  +  {DET}T[M3000]{RE} 

6  I- 

+  (ETAKJ  }TCM400K0]{UI  }{SIJ) 

-  { ET AKJ }T ( ( { ET ALI }T [M500KOL] ) { KDI FI J } ) {RK} 

+  v|j (ETAKJ }T[M40KL0]{UBARLKRK}  j+1j  (62) 

(FRE)e  =  (DET  }T[M3000]  +  (ETAKJ  }T[M40K00](DET}J{RE}j+1 

+  ^  -{UBARK}T[M30K0](RE}  +  {ESK}T[M40000](DET}{RE} 

+  {ETAKJ}1  [M400KO](UI  XSIJ'ESK} 

-  (ETAKJ }T( ( (ETAKI }T [M500 KOL } ) { ED I  FI J > ) {RE} 

+  { ETAKJ }T[M40KL0](UBARL}(RETJj+1 j  (63) 
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A  few  comments  regarding  notational  structure  are  appropriate.  The 
elemental  hypermatrices  [M...],  for  an  n-dimensional  problem  description, 
are  evaluated  once  and  for  all  by  integrals  over  a  master  finite  element 
domain  R^.  For  a  one-dimensional  domain  R1,  i.e.,  [A...],  all  matrices 

through  [A4 _ ]  are  listed  in  Appendix  A,  for  both  the  linear  (k=l)  and 

quadratic  (k=2)  cardinal  basis  {N^}.  For  two-  and  three-dimensional  domains 
Rn,  ie.,  [ B . . . ]  and  [C...],  all  element  matrices  through  [B4...]  and  [C2...] 
are  listed  in  Appendix  B  and  C,  respectively,  for  the  bi-linear  tensor 
product  cardinal  basis  (Nil.  The  discrete  indices  J,  K,  L,  occurring  in 
both  matrix  and  variable  (FORTRAN)  names,  are  tensor  summation  indices  with 
range  1  <  (J,K,L)  <  n.  The  multi-pole  coefficient  ^2  is  expressed  in  terms 
of  Cartesian  scalar  compoennts  v^j,  with  distinct  values  for  each  dependent 
vari  able. 


The  arrays  (DET }  and  {ETAKJ}  contain  nodal  values  of  the  determinant 
J  and  elements  of  Sn^/SXj.  The  nodal  values  of  the  contravari ant  scalar 
components  of  the  convection  velocity  u^,  equation  39,  are  denoted  {UBARK}. 

The  elements  of  {SIGIJ}  are  nodal  values  of  the  stress  tensor  computed  in 
principal  coordinates  in  terms  of  u  ■  =  m-/p.  Equation  60  is  evaluated  only 

J  v 

for  the  laminar  flow  Stokes  stress  tensor  definition,  equation  4,  and  Sjj 
is  the  Kronecker  delta.  The  additional  terms  resulting  from  inclusion  of 
the  Reynolds  stress  constitutive  equation  9,  see  equation  8,  are  readily 
evaluated  using  the  illustrated  expansions.  In  equation  61,  (GSR)  denotes 
nodal  values  of  e  5  g/p,  while  (UIUI }  contains  nodal  values  of  the  specific 
kinetic  energy  £u-u..  For  equations  62-63,  {KDIFIJ}  and  { ED  I  FI J }  contain 
as  elements  the  nodal  values  of  the  "effective  diffusion"  coefficient  tensors. 


Ckpuiuj  k/e  -  pS.jj  and  C£pu|uf  k/e,  respectively.  Further  in  equation  63, 
the  elements  of  {ESK}  and  {SIJ-ESK}  are  nodal  values  of  tie/Y,  and  C*pu "u^'e/k. 

D  b  b  J 

In  all  equations,  {*Kj+1  denotes  -  (* )j •  The  notation  ]j+1  j  defines 

evaluation  of  the  argument  at  t j  ,  and  at  t j ,  followed  by  addition  after 
multiplication  by  0  and  (1-0),  see  equation  24. 


5.  Tensor  Jacobian  Matrix  Statement 

The  algorithm  matrix  statement  is  completed  by  construction  of  the 
tensor  matrix  product  Jacobian,  see  equations  27  and  42.  Equations  56-63 
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provide  the  expressions  for  { FI >e ,  which  are  analytically  differentiable 
by  members  of  {QJ >e .  The  construction  of  certain  terms  involves  differ¬ 
entiation  with  respect  to  the  parameter  =  rn^/p.  Using  the  chain  rule, 
and  equations  38-39, 


9  _  3  ,3  3{MK}  ,  3  3 { U K } 

mue  "  JWT  bTrkT  HmTT  MukT  sTmTF 


3  _ 

3WI 


+  det  J 


3  1  3 

1W  +  P  3iDKi 


(64) 


3  _  3 

f\ 

3 

3lRJe  3{R)  ' 

3{Qk> 

(65) 


Denote  as  the  element  average  of  det  J  (Sr^/SX-j  )e ,  and  define  K  to 
signify  the  discrete  free  index,  ie.,  no  summation  implied,  corresponding 
to  3/3r|,,  see  equations  53-54.  The  non-empty  tensor  product  Jacobians  for 
equations  56-61,  suppressing  the  subscript  e  throughout,  are 


[JRR]e  =  {DETK}T[A3000]  +  vj^{  ETAKJ  }T[A40K00]{DETK} 

T  fm 

{ETAKJ}  [M40KL0]{UBARL}  - 


.At  2 
+  T  VRJ 


P2J 


{ETAKJ }T[M40K0L]{R} 


[JRMIL  = 


At 


•{ETAKI}T[A30K0]  +  v^Dki{ETAKJ}T[A30KK] 


(66) 


[JMIR]e  =  f- 


At 

2 

-2 

P' 

{MI}T[A30K0]  +  vf j (MI }TCA4KK00] {ETAKJ } 


[JMIMI]e  =  {DETK}T[ A3000]  +  vjj{ETAKJ}T[A40K00]{DETK} 


At 


■{UBARK}T[A30K0]  -  Dki {MI }T[A30K0]  i 

iTr 


+  v? t ( {UBARK} 1  [ A40KK0] {ETAKJ}  +  4-{MI}‘  [A4KK00] { ETAKJ } ) 

lu  P 

At 


[JMIMJ]e  -  ^  Dkj  |4 


■{MJ }T[ A30K0]  +  vj  L 01 J >T  C  A4KK00 ] { ETAKL } 


[JMIP]e  =  -  ^r{ETAKI }T[A30K0] 
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[JMISIJ ]e 

[JGft]e 

[JGMI]e 

[OGG]e 

[JGP]e 

[JGSIJ]e 

[JGQJ]e 

[JpR]e 

[JPMK]g 

[JPG]e 

[GPP]g 

[JSIJR]e 


=  ^{ETAKJ}T[A40K00]{DETK} 


(67) 


At 

2 


\ 

W 


{G+P  }T[  A30K0]  +  vg  j  (  G  }T  [  A4KK00  ]  {  ET  AKJ  }[ 


J 


At 


2  Mi1 


■  {G+P}T[A3OK0]  +  vGL{G}T[A4KK00]{ETAKJ>j 
=  {DETK}T[A3000]  +  Vgj{ETAKJ}T[A4OK00]{DETK} 
AtQlJBARK}T[A3OK0]  +  v^{UBARK}T[A40KKO]{ETAKJ} 


4^{UBARK}T[A30K0] 


=  -y-{ETAKJ}T[A40K00]{UI} 
=  ^  Dlj{ETAKL}T[A30KO] 


(68) 


full 

pkl 

l  2 

H 

'ill' 

2 


{DETK}7 ( [A40000]{MK> 


{DETK}T([A40000]{UK})  +  fi|{DETK}T([A40000]{M!<}) 


=  - (y-1 ) {DETK}  [A3000] 


{DETK}  [A3000] 


(69) 


_  y 


m.  {ET  AKJ  }T[  A30K0]  +  m  .{ETAKI  }T[A30K0] 

... '  J 


_  i 

3 


6I(Jm£{ETAKL}T[A30K0] 
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[JSIJMI]e  =  -  #  {ETAKJ}T[A30K0]  -  f6IJ{ETAKL}T[A30K0] 

[JSIJSIJ]e  =  {DETK}T[A3000]  (70) 

-  -2 

r—  j  m.  T  — 

[JQIR]  =  £  g{ETAKI}  [A30K0]  -  —{ETAKI } 1  [A30K0] 

Krn  •  j 

[JQIMI]  =  -=j-  {ETAKI }'[A30K0] 

C  P 

[JQIG]  =  -  -i-{  ETAKI  }T[A30K0] 

G  p 

CJQIQI]e  =  {DET  K}T [A3000]  (71) 


The  scalars  ,  p  and  g  are  element  average  values  of  {MI},  {R>  and 
{G>  respectively.  In  the  formation  of  equations  66-71,  the  Boolean  index 
K(=l)  in  various  A-matrices  has  been  permuted  to  facilitate  differentiation 
of  each  expression  by  the  last  right  contraction  matrix.  There  is  consid¬ 
erable  commonality  pervading  the  individual  element  Jacobian  constructions. 
The  elements  of  each  Jacobian  are  formed  on  the  finite  element  domain  R1, 

G 

using  the  element  matrices  listed  in  Appendix  A,  and  then  assembled  into  the 


global  form  using  the  operator  Sg,  see  equation  55.  In  actual  practice,  the 
column  matrix  {<5QI}  is  ordered  on  degrees  of  freedom  at  a  node,  e.g.. 


T 

{...,  6Rj,  6Uj,  6Ej,  <5Pj ,  <$Rj+1,  •••}  •  Hence,  the  global  tensor  Jacobian 
[Jri],  see  equation  42,  is  a-block  tri -diagonal .  using  the  linear  (k=l)  finite 


element  basis,  and  a-block  penta-diagonal  for  the  quadratic  (k=2)  cardinal 


basis. 
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6.  Comments  on  the  Algorithm 


Equations  56-63  are  an  exact  statement  of  the  elemental  matrix  equiva¬ 
lent  of  the  finite  element  algorithm  statement,  equation  24,  which  upon 
assembly  on  UR^  yields  equation  24.  While  the  master  hypermatrix  symbology 
[MP...]  is  independent  of  the  completeness  of  the  approximation  subspace, 
ie.,  the  degree  k  of  {^(n)},  equation  21,  the  order  of  each  square  matrix 
is  equal. to  (k+l)n.  In  addition,  since  each  master  matrix  is  a  hypermatrix 

D 

of  degree  one  at  least,  P  >  1,  there  are  (k+1)  entries  for  each  element  of 
each  matrix  [MP. . .]  . 

The  associated  inner  product  DO-loop  to  form  {FI},  especially  for 
k>  1,  would  be  of  excessive  length  on  a  scalar  machine.  The  algorithm  matrix 
statement  can  be  simplified  in  this  regard,  with  the  commission  of  interpola¬ 
tion  error  only.  Each  Boolean  index  "0"  appearing  in  a  hypermatrix  [MP...] 
signifies  the  matrix  inner  product  with  the  element-dependent  column  (row) 
matrix  {•}  corresponding  to  an  interpolation  on  Rn.  For  example,  the  lead 
term  in  each  of  equations  56-63  is 

{DET}g[M3000]e 

In  the  instance  of  the  linear  transformation  x  =  n»  The  elements  of  {DET} 

T  e 

are  each  equal  to  a  constant  times  the  element  measure.  Hence,  {DET}  = 

T  ® 

(DET)e{0NE}  ,  where  the  elements  of  {ONE}  are  unity,  which  yields 

{DET }g[M3000]  =  (DET)e{0NE}T[M3000] 

=  (DET)e[M200]  (72) 

where  (DET)e  is  the  scalar  measure  (area,  volume)  of  the  finite  element  domain 
Rg.  As  presented  in  Section  V,  for  k=l,  n=2,  and  a  general  quadrilateral 
domain  R^,  the  elements  of  {DET}g  are, 
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®ET}e  =  } 


e 


(73) 


f (X2-X1) (Y4-Y1) 
(X2-X1) (Y3-Y2) 
(X3-X4) (Y3-Y2 ) 
[  (X3-X4 ) (Y4-Y1) 


(X4-X1) (Y2-Y1) 
(X3-X2 ) (Y2- Y1 ) 
(X3-X2 ) (Y3-Y4) 
(X4-X1) (Y3-Y4) 


with  counterclockwise  node  numbering  starting  in  the  lower  left  corner  of 
Rg,  see  Figure  la).  Defining  the  sum  of  the  entries  as  (DET)g,  and  for  a 
"decent"  element  aspect  ratio,  equation  72  can  be  approximated  as 


(DET}g[B3000]  5  (DET)e[B200] 


(74) 


with  the  commission  of  interpolation  error  only.  This  operation  reduces  by 
(k+1)^  the  number  of  multiplications  to  form  this  term  in  {FI}. 

In  this  context  then,  every  hypermatrix  inner-product  statement  in 
equations  56-63  is  eligible  for  reduction  of  hypermatrix  degree  through  ele¬ 
mental  averaging.  In  terms  of  error,  thi s  i s  appropriate  only  for  element- 
dependent  data  that  are  sufficiently  smooth.  In  general,  this  rules  out  the 
averaging  of  any  dependent  variable,  since  non-smooth  solutions  must  be  ad¬ 
mitted.  However,  for  "decent"  discretizations  UR^,  the  geometric  data  should 
be  eligible  for  averaging  on  sufficiently  refined  grids.  Hence,  for  example 
in  equation  57, 

{ ET  AKJ  }T[M40K00  ]  {DET}_  ■+■  (ETAKJ)q(15ET).[M2KO] 

G  G  G 

{ ET  AKJ  }  J  [M40  KLO  ]  {UBARL  }Q  +  (ETAKJ)  Q{UBARL}T[M30I<L]  (75) 

G  6  G  G 

In  equation  62,  for  example 

{ETAKJ )I ( ( (ETALI }J[M500K0L3) {KDIFIJ}  ) 

G  G  G 

(EFAKJ)  (ETAKL)  {KDI  FI  J } J  [M30  KL]  (76) 

c  G  G 


Hence,  no  hypermatrix  of  degree  P  >  1  is  required  formed  or  stored,  provided 
the  discretization  of  R^  is  of  sufficient  quality.  The  averaging  of  the 
grid  data  in  the  Jacobian  matrix  formulation,  equations  66-71,  is  thus  also 
appropriate. 
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A  second  point  to  note  is  that  definition  and  use  of  the  stress  tensor 

a..,  and  heat  flux  vector  q.,  as  dependent  variables  has  permitted  algorithm 
"I  J  J 

formulation  without  evaluation  of  second  order  spatial  derivatives  in  the 
generalized  coordinates  framework.  Hence,  the  troublesome,  sometimes 
destabilizing  mixed  derivatives  resulting  from  viscosity  terms  in  the  direct 
(laminar  flow)  formulation,  are  totally  absent.  The  penalty  for  the  depen¬ 
dent  variable  construction  is  a  significant  increase  in  size  of  the  block- 
banded  tensor  product  Jacobian.  Conversely,  a  general  turbulent  flow  predic¬ 
tion  can  ostensibly  be  handled  with  ease,  by  expansion  of  the  defining  equation 
for  a.,  and  q..  The  fact  that  the  algebraic  "constitutive"  equations  for  p, 
a-  •  and  q.  are  handled  directly  within  the  weighted  residuals  algorithm  lends 

'  J  J 

an  overall  uniformity  that  simplifies  construction. 


Finally,  non-iterative  and  direct  steady-state  forms  of  the  algorithm 
are  special  cases  of  the  formulation.  The  non-iterative  construction  simply 
constitutes  acceptance  of  the  first  solution  (SQI}j+^  of  the  Newton  algorithm 
using  only  evaluation  of  (FI > ,  see  equation  25.  Since  the  iteration  index 
p  >  0  by  definition,  (QI)j+^  =  {QI }j .  Therefore,  (QI}J+1  =  (0}  in  equations 

56-63,  and  the  corresponding  expressions  in  brackets  are  not  evaluated. 

Furthermore,  At/2  ->•  At  and  the  evaluation  •].  ,  .  reduces  to  •]..  Since 

J  J*  >  J  J 

the  terms  in  (FI)  involving  v* .  have  been  eliminated,  so  are  the  corresponding 
terms  in  the  tensor  product  Jacobians,  equations  66-68.  The  multiplier 
At/2  remains  appropriate  in  [Jr,],  and  evaluations  are  made  using  = 

{QI}..  This  procedure  thus  reduces  the  algorithm  operations  count  by  a 

vJ 

significant  factor,  at  the  expense  of  removal  of  v1 •  from  the  construction. 

aJ 

The  direct  steady-state  algorithm  is  identical  to  the  non-iterative 
formulation  except  that  At/2  -*  At  in  the  Jacobians.  This  multiplier,  common 
to  all  elements  of  {FI},  can  be  divided  out,  yieTdingAt  1  as  a  scalarmultiplier  on 
{DET }T[A3000]  in  the  self  coupling  Jacobians  [JQQ]  for  the  initial-value  depen¬ 
dent  variables.  This  yields  the  tensor  matrix  generalization  of  the  "approxi¬ 
mate  factorization"  procedures  devised  for  finite  difference  algorithms,  cf. 
[5,6]. 
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SECTION  IV 


THEORETICAL  ANALYSIS,  ACCURACY  AND  CONVERGENCE 
1 .  General  Concepts 

The  theoretical  analysis  of  accuracy  and  convergence  aspects  of  a 
finite  element  algorithm  is  typically  quantized  in  inequalities  in  Sobolev 
norms.  The  Navier-Stokes  problem  class  for  finite  Reynolds  number  is  non¬ 
linear  elliptic  with  initial-value  character.  For  infinite  Reynolds  number 
(inviscid  flow),  the  resulting  (Euler)  equation  set  is  non-linear  hyperbolic. 

The  theoretical  analyses  are  exact  only  for  the  linear,  one-dimensional  initial- 
value  problem  description,  and  the  linear  elliptic  multi-dimensional  problem 
class,  cf.  [11].  The  convergence  statements  are  constructed  by  quantizing 
the  total  approximation  error  in  terms  of  component  parts  as: 

semi-discrete  approximation  error:  c^(x,nAt)  =  q(x,nAt)  -  qh(x,nAt) 

discrete  approximation  error:  a(x,nAt)  =  q(x,nAt)  -  Q(x,nAt) 

temporal  truncation  error:  t(x,nAt)  =  Q(x,nAt)  -  qh(x,nAt)  (77) 

For  any  choice  of  norm,  the  triangle  inequality  yields 

|Ie||  =  || a  +  x ||  <  ||al|  +  ||t||  (78) 

The  semi-discrete  approximation  error  in  a  finite  element  solution  of 
a  linear  parabolic  equation,  using  the  fully  explicit  Euler  integration 
algorithm,  0=0  in  equation  24,  satisfies  the  inequality  [11,  p.  408]. 

|| e  (nAt) |||^Rij  <  CiAe^  ^ ||q (nAt)  ||^k+l ^  +  CzAt ||Q0 ^  (79) 

The  Ca  are  constants  independent  of  Ae,  the  measure  of  the  uniform  finite  ele¬ 
ment  discretization  of  R1,  m=l  for  the  parabolic  specification,  and  ||Q0||^m 
is  the  tf1  norm  of  the  interpolation  of  the  initial  data. 
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The  corresponding  statement  for  smooth  solutions  of  a  linear  hyper¬ 
bolic  equation,  and  using  explicit  Euler  integration,  is  [11,  p.  418] 


h 


(nAt)||  <  C1Ag+1||q(nAt)||Hk+l^Rl^+ C2At||Q0||^o  +  C3AeJ^  ||q(T)||dx+i  (80) 


nAt 


The  remain  constants  independent  of  Ae,  the  first  term  is  the  familiar 
bound  on  discretization  refinement,  and  the  second  term  is  the  iH°  norm  of 
the  initial-data  interpolation.  The  third  term  is  an  additional  discretiza- 
tion-sensi ti ve  error  generated  through  propagation  of  the  k+1 —  derivative 
of  the  exact  solution  over  nAt.  The  exponent  of  unity  on  Ag  is  always 
smaller  than  the  exponent  k+1  on  the  lead  term.  Therefore,  enhanced  conver¬ 
gence  through  use  of  a  more  complete  basis,  k  >  1  in  equation  21  is  not 
expected. 

The  convergence  proof  for  a  linear  elliptic  boundary  value  problem 
definition,  on  a  multi -dimensional  domain,  can  be  written  as  [11] 


II e  Hh06 (Rn )  -  CiaY Il9 Hhp (3R) 

where  ||g||^p  is  the  Sobolev  norm  of  the  boundary  data  specification  of 
the  problem.  Cx  remains  a  constant  independent  of  the  extremum  mesh  measure 
A,  a  <  1 ,  and  y  =  min[k,  p  +  s  -  4,  1  -  a]  for  p  >  4  and  s  >  4.  The  parameters 
p  and  s  pertain  to  the  regularity  of  the  discretization  of  UR",  and  conver¬ 
gence  can  become  independent  of  the  degree  k  of  {N^}  with  distorted  meshes 
having  large  discretization  aspect  ratios.  Otherwise,  y  takes  on  the  value  k. 

In  regards  to  accuracy,  an  important  additional  theoretical  aspect 
is  extremization  of  the  (semi-)  discrete  approximation  error  by  a  finite 
element  solution.  Strang  and  Fix  [24]  prove,  for  the  one-dimensional 
Laplacian,  that  the  corresponding  linear  (k=l)  finite  element  solution 

h  h 

produces  an  absolute  minimum  discrete  approximation  error  e  =  q  -  q  , 

l_  L. 

measured  iin  the  energy  norm  E(e  ,en),  which  is  defined  as 
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E(qh-qh)E.fRi^<lx  (82) 

The  theorem  states, 

E(q-qh,  q-qh)  <  E(q-q»  q-q)  (83) 

for  all  possible  choices  of  an  alternative  approximation  q,  including  the 
(linear)  interpolation  of  the  exact  solution  and/or  a  finite  difference 
solution.  Generalizing  to  degree  k,  and  using  orthogonal ity,  they  further 
prove  that  the  energy  in  the  approximation  error  is  bounded  in  the  form 

E(e\eh)  S  C,A2k||q||2M  (84) 

Again,  Cj  isa constant  independent  of  the  mesh,  A  is  the  extremum  measure, 
and  ||q ||2k+l  is  the  square  of  the  k+1  Sobolev  norm  of  the  exact  solution. 

~  2  k+1  ^ 

Hq|lHk+1  -  |Rl  qZ  +  (a?)  +  +  (i?+ij  dx  (85) 

The  statement  corresponding  to  equation  81  is  then, 

E(e\Sh)  <  C2A2k||g||2H„  (86) 

where  ||g||^o  is  the  H°  Sobolev  norm  of  the  boundary  data,  which  need  only 
be  square  integrablc. 

In  approaching  the  Navier-Stokes  problem  class.  Baker  and  co-workers 
document  results  of  detailed  computational  experiments,  defined  to  confirm 
the  convergence  exponent  2k,  equation  86,  or  k  in  the  equation  79,  for  select 
model  equation  systems.  For  a  linear,  transient  parabolic  heat  conduction 
equation  [25],  the  convergence  data  were  exactly  interpolated  by  k  over  the 
range  1  <  k  <  3.  This  confirmation  was  similarly  obtained  [26]  for  1  <  k  <  2, 
for  solution  of  the  mildly  non-linear  laminar  flow  boundary  layer  equation 
system,  fora  range  of  pressure  gradients.  Significantly,  the  energy  in  the 
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h  h 

semi -discrete,  k=l  approximation  error  E(e  ,en),  equation  86,  was 
confirmed  extremum  in  comparison  to  the  Crank-Nicol son  algorithm  solution, 
including  use  of  non-uniform  discreti.'2'ations.  Confirming  results,  on  error 
extremization  and  convergence  rate  exponent  for  1  <  k  <  2,  are  reported  [27] 
for  the  highly  non-linear  turbulent  boundary  layer  equation  system  solution. 

Of  principal  additional  significance,  the  energy  norm  was  a  strongly  nonlinear 
function  of  the  turbulent  "eddy  viscosity,"  yet  the  semi-discrete  error 
convergence  data  were  closely  interpolated  by  lines  of  constant  slope  equal 
to  2k.  Finally,  for  a  smooth  solution  to  a  linear  hyperbolic  equation  [25], 
the  convergence  rate  was  indeed  independent  of  k,  confirming  dominance  of  the 
third  term  in  equation  80.  In  all  cited  cases,  on  a  sufficiently  refined 
discretization,  the  absolute  level  of  semi-discrete  approximation  error  for 
the  k  >  1  algorithm  solutions  was  uniformly  less  than  that  for  the  k=l  solu¬ 
tions.  Furthermore,  in  comparison  to  the  equal  complexity  Crank-Nicol son 
finite  difference  algorithm,  the  k=l  finite  element  semi-discrete  approxima¬ 
tion  error  in  ||-  ||^i  was  uniformly  mimimum  for  any  discretization  of  R1. 

2.  The  Parameter  Set  1i? 

A  principal  requirement  for  high  Mach  number  Navier-Stokes  or  Euler 
solutions  is  prediction  of  non-smooth  solutions,  ie.,  ones  exhibiting  a 
finite  number  of  finite  discontinuities  (shocks).  The  finite  element  theo¬ 
retical  literature  does  not  contain  convergence  statements  analygous  to  equa¬ 
tions  79-80,  nor  is  extremization  of  the  semi-discrete  approximation  error 
guaranteed.  The  theoretical  requirement  is  to  estimate  a  reasonable  value 
for  the  multi-pole  expansion  coefficient  set  ~$2>  equation  22.  For  the  one¬ 
dimensional  form  of  the  continuity  equation  1,  and  for  an  imposed  constant 
convection  velocity  u-  =  U0 ,  the  exact  Fourier  solution  is, 

J 

q(x,t)  =  Q0  exp[ia)(x-U0t)]  (87) 

corresponding  to  propagation  of  the  initial  data  Qq  with  group  velocity  U 
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The  semi '-discrete  Fourier  solution  is, 


qh(x,t)  =  Qq  exp  [iw(jAx  -  rUQt)] 


(88) 


where  r  =  a  +  is,  where  a  and  6  are  real  numbers,  i  =  v^I,  w  =  2-rr/A  is  the 
wave  number  for  Fourier  mode  of  wavelength  A,  and  x  =  jAx,  j  =  0,1,2,..., 
is  represented  by  discrete  intervals  of  (uniform)  measure  Ax.  Raymond  and 
Garder  [28],  for  the  definition  $2  =  vAxi,  where  v  >  0  is  a  scalar  parameter, 
determined  the  k=l  algorithm  functional  forms  for  a  and  6  as. 


a 


1  + 


+  0 (ds ) 


(89) 


' 6  =  “  T?  d3  +  °(dS) 


(90) 


In  equations  89-90,  d  =  ooAx  =  27T/n,  where  n  is  the  discrete  Fourier  mode 
index,  Af,  =  nAx,  and  0  indicates  order.  Hence,  the  semi-discrete  Fourier 
solution  q*1  can  be  made  a  sixth-order  accurate  approximation  to  equation  87 
by  eliminating  the  0(d4)  term  in  equation  89,  yielding  v  =  (15)~2.  Corre¬ 
spondingly,  6  <  0  in  equation  90,  and  an  artificial  dissipation  mechanism 
becomes  introduced,  yielding  the  semi-discrete  approximation  as 

qh(x,t)  =  Q0exp[iw(jAx-Uot)]exp[-w4kt  +  0(Ax6)]  (91) 

Here,  k  s  U  (Ax^l&TF  is  the  coefficient  of  "artificial  viscosity,"  and  the 
dissipation  mechanism  is  highly  phase  selective  due  to  the  to4  factor. 

The  original  analysis  was  expanded  [15]  by  redefining  the  dissipation 
parameter  §2  in  the  form, 

tz  =  Ax[v16t  +  v26x]i  (91) 

where  <5.  and  6  are  Kronecker  delta-type  functions,  yielding  v1  operating 

L  A 

on  the  time  derivative,  and  v2  operating  on  the  spatial  derivative  term,  only, 
in  the  one-dimensional  form  of  equation  1.  Proceeding  through  the  substitu¬ 
tions  yields,  for  the  k=l  algorithm. 
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a  =  1  -  d2(v1  -  v2 )  +  d4  -  Y^n  +  +  (v1  -  v2)  (v1) 


v 


6  =  d 

v1  -  v2 

-  d3 

_  _ 

L 180  12 

gr-  (v1  -  v2)(v1)' 


+  0(d6) 


+  0(d5) 


For  the  quadratic  (k=2)  algorithm,  the  form  of  equation  92  is  [15]. 


a=l-  4v1v2  +  d2 


+  dl 


-14  +  184v1v2  -  SOCv1)2  +  240v2(v1) 


/15 


2  5 42  i  2  i  8/  i \ ^  672  2/  i\^ 

T-T50vv  +3(v)  -TTv(v) 


+  16 (v1)  -  64v2(vx) 


+  0  (d6 ) 


(92) 

(93) 


(94) 


Setting  v1  =  v  =  v2  in  equations  .92-93  yields  the  results  of  the 
original  analysis  [28].  Enforcing  sixth-order  accuracy  for  the  k=l  algorithm 
construction,  equation  92,  yields  the  constraint 


d2 

180  '  <vl>* 

’  +  (v1)2 

d2 

|  -  H 

+  V1 

(95) 


Figure  2  is  a  plot  of  equation  95  with  n  as  a  parameter.  Sixth-order  accuracy 
can  be  achieved  only  for  v1  >  0;  for  any  level,  v2  ranges  over  an  order  of 
magnitude  dependent  upon  n,  with  the  largest  levels  required  for  the  shortest 

_JL 

wave  lengths.  All  curves  converge  to  the  point  v1  =  (15)  2  =  v2.  Defining 
v1  =  0  renders  both  the  k=l  and  k=2  algorithm  construction  for  the  semi-dis¬ 
crete  approximation  q*\  equation  91,  a  second-order  accurate  representation  of 
the  analytical  solution. 

The  phase  velocity!  of  each  Fourier  mode  can  be  determined  for  the 
fully  discrete  solution  Q,  defined  as 


Q  =  q]?  (Ax ,  nAt )  =  gnexp[iwjAx] 


(96) 
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Substituting  equation  96  into  the  one-dimensional  continuity  equation,  inserted 
into  the  algorithm  statement,  equation  22,  integrating  using  the  implicit 
trapezoidal  rule,  0  =  i  in  equation  24,  and  defining  B2  =  vAx? ,  yields  the 
amplification  factor  in  equation  96  in  the  form  [29]. 

1+4  cos(caAx)  -  3Cv  sin2(ia)Ax)  -  i-|(C+2v)sin(wAx) 

9  ~  -5 

1+4  cos(wAx)  +  3Cv  sin2(4u)Ax)  +  i-|-(C=2v)sin(wAx)  v  1 

In  equation  97,  C  e  U0Ax/At  is  the  non-dimensional  time  step  (Courant  number), 
and  v  is  the  dissipation  parameter.  For  v  =  0,  |g|  =1  for  all  At  and  Ax, 
indicating  the  k=l  finite  element  algorithm  will  neither  damp  nor  amplify 
solution  disturbances,  ie.,  it  is  neutrally  stable.  For  all  v  >  0,  |g|  <1 
and  an  artificial  dissipation  mechanism  becomes  introduced. 

The  non-dimensional  phase  velocity  0^  of  each  discrete  Fourier  mode 
is  defined  in  terms  of  the  real  and  imaginary  components  of  g  as  [29] 


0 


Imag 

& 

~ 

real 

!g! 

(98) 


Figure  3  graphs  equation  98  for  select  levels  of  C  and  v,  as  a  function  of  the 
discrete  wavelength  parameter  n,  recall  Xn  =  nAx.  The  exact  solution  is 
0n  =  1  for  all  n.  The  fully  discrete  solution  cannot  resolve  the  "2Ax"  mode,  n=2. 
For  C  <  0.5  and  v  =  0,  0n  >  0.95  for  all  n  >  5,  ie.,  the  discrete  solution 
will  artificially  disperse  these  low-wave  length  solution  modes  yielding  solu¬ 
tion  "wiggles,"  cf.  [17].  Increasing  the  Courant  number  further  retards  the 
phase  velocities  for  all  n.  Setting  v  =  (15)~*  and  C  =  0.5,  yields  On  >  0,97 
for  all  n  >  2,  which  essentially  eliminates  the  short  wave  length  phase 
dispersion.  Unfortunately,  this  excellent  phase  accuracy  has  associated  with 
it  a  high  level  of  artificial  dissipation  which  excessively  diffuses  the 
gradients  in  the  physical  solution.  For  comparison,  the  dashed  curves  in 
Figure  3  interpolate  the  data  for  ©n  as  obtained  for  the  Crank-Micolson  second 
order  accurate,  0  =  4  finite  difference  solution  algorithm  applied  to  the 
one-dimensional  continuity  equation  with  UQ.  Unacceptable  phase  dispersion  is 
indicated  for  all  n  <  20. 
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PHASE  VELOCITY  0 


The  multi-dimensional,  multi-dependent  variable,  generalized  coor¬ 
dinates  definition  for  the  dissipation  parameter  3“  is 

"$2  =  (det  J)[ P  6f  +  v2  6  ]  (99) 

CC  U  X 

In  equation  99,  det  J  is  the  measure  of  a”,  subscript  a  denotes  the  appro¬ 
priate  member  of  qj,  and  the  parameter  vectors  vK  are  expressed  in  terms  of 

cc  ct 

scalar  components  v£.  in  the  x .  coordinate  system,  see  equations  56-63. 


Figure  3.  Fourier  Phase  Velocity  Distribution  for  Finite  Element  and 
Finite  Difference  Algorithms. 
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SECTION  V. 


RESULTS  AND  DISCUSSION 
1.  Quasi-One  Dimensional  Flow 

Several  essential  aspects  of  algorithm  performance  regarding 
accuracy,  convergence,  etc.,  can  be  quantized  by  examination  of  quasi- 
one-dimensional  inviscid  flows.  Defining  the  convection  velocity  u  5  m/p, 
and  the  flow  cross-sectional  area  as  A(x),  the  governing  differential  equa¬ 
tion  set  is  formed  from  equations  1-7  as 


L(p) 

L(m) 

L(g) 

L(p) 


■  3t  +  |xw +  3S"A  ■ 0 

(100) 

■  Tt  *  lx [um  +  + «  *£* 

=  0 

(101) 

=  It  +  lx  ^ug  +  upl  +  u(g+p^ 

3£nA  _  q 

3x 

(102) 

=  p  -  (y-l)Cg  -  ium]  =  0 

(103) 

The  finite  element  algorithm  statement,  equations  56-59,  requires 
the  metric  data  {DETl  and  {ETAKJ}  .  For  the  affine  coordinate  transfor- 

G  G 

mation  x1  =  ru,  and  with  the  origin  of  the  oi  coordinate  at  the  element 
centroid,  cf.  Figure  1,  the  elements  of  the  cardinal  bases  {N^rOK 
1  <  k  <  2,  are 


(104) 


{N2 (n ) )  =  i 


-n(l-n)  ] 

(l-n)(l+n) / 

n(l+n)  J 


(105) 
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By  definition,  Xi  =  {N^ (n) }"^{XI  }g ,  and  det  J  =  det[3xi/3m]  •  Denoting 
the  elements  of  {XIl  as  the  left  (L),  right  (R)  and  middle  (M)  x.-coor- 
dinates,  where  M  =  i(R+L)>  then 

k=l:  DETe  =  i(R-L) 

k=2 :  DETe  =  i(R-L)  -  -^(R+L)  (106) 


Defining  the  element  measure  as  A  =  R-L,  then  DET  =  A  „/2  for  both  the 
linear  and  quadratic  basis.  By  the  same  procedures,  the  sole  non-vanishing 
element  of  {ETAKJ>e  is  in  the  (1,1)  location  with  a  value  of  unity. 

Therefore,  the  finite  element  algorithm  statement  becomes  quite 
simplified  for  the  quasi -one-dimensional  situation.  The  specific  form  for 
equations  56-59  is 


{FR)e  =  |  *e[A200]  +  vi[A210  Jj 


{R>j+l 


,  At 
+  -T 


[A2 10] {M>  +  vi{U) ' [A301 1 ] (R) 


+  (A)  ( [A41000] {U> ) {R} 


j+l»j 


;  =  {0} 


(107) 


{FM}e  = 


A  [A200]  +  vi[A210]) 


At 


{V}T[A3010]{M}  +  [A2 10 ] {P }  +  v I { U }T [A3011]{M> 


+  {A}  ([A41000]{U}){M} 


j+1  J 


(0) 


(108) 


(FG}e  - 


Ae  [A200]  +  v13[A210] 


{G}' 


+1 


At 


{ U}T [A30 10 ]  [{G>  +  {P}|  +  { U}T [ A301 1  ] { G } 
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+  {A}T ( [A41000] {U} ) ({G}+{P}) 


j+l>j 


•  =  {0}  (109) 


(FPL  =  A o[A200]{P}  -  A.(y-l) 

e  c  c 


[ A200] {G}  -  |  {U}T;[A3000]{M} 


=  {0} (110) 


The  Jacobian  contributions  for  the  Newton  algorithm  are,  recall  equations 
66-69. 


8(FR} 

wt~ 


At 


[JRR]  =  Ae[A200]  +  vJ[A210]  +  ^ 

-  -~{R}T[A3110]]  +  { A}T[ A4 1000] 


v?  {U}T [ A301 1 ] 


{U}-  rr(R} 


[JRM]  =  [A210]  +  Vi[A211]  +  {A)T[A3100] 


) 


j+1 


[JRG]  =  0  =  [JRP] 


(111) 


3  (  FM  } 
aTRT 


=  [JMR]  = 


[A3010]  +  v*[A3110] 


+  [A40001] {A} 


P 

j+1 


[JMM]  =  Ao[A200]  +  v*[A210] 


^{U}T|Ta3010]  +  [A3011]  +  [A4000 1  ] { A}  P 


At 

~T 


j+1 

-TP 


-=•]  {M}T[[A3010]  +  vf[A3110]  +  [A40001]  {A}] 

lp)  L  Jj+1 


[JMG]  =  [0] 


[JMP]  =  4f{A210] 


(112) 


=  [JGR]  =  -  ~||~J|j[G}T([A3010]  +  Vg[A3110]  +  [A40001]{A}) 

T  ~ rP 

+  {P }  [A3010] 

Jj+1 
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[JGH]  -  f  (i) 


{G}  ( [A30 10]  +  v$[A3110]  +  [A40001]{A}) 


+  {P}  [A3010] 

-Jj+1 


[JGG]  =  A_[A200]  +  v*[A210] 

t  * 


•+  fmJ 


[A3010]  +  v|[A3011]  +  [A40001] {A} 


-tP 


J+l 


[jgp]  =  fmJ 


[A3010]  +  [A40001] {A} 


P 

j+l 


(113) 


_1_  9{FP} 
Ag  3{R} 


E  [JPR]  =  - 
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rVl (M}T[A3000] 


P2J: 


P 

j+l 


[JPM]  -  -^2“ 


|{li)}T[A3000]  +  (~|{M}T[A3000]lP 


~^j+l 


[JPG]  =  - (y~ 1 ) [A200] 


[JPP]  =  [A200] 


(114) 


In  equations  112-114,  the  superscript  bar  on  m  and  p  indicate  the  element  i 
average  values.  The  elements  of  (A>e  are  the  nodal  values  of  £n  A(x).  The 
defined  standard  matrices  and  hypermatrices  are  listed  in  Appendix  A  for 
both  cardinal  basis  formulations,  1  <  k  <  2. 


2.  Accuracy  and  Convergence,  Riemann  Shock  Tube 

The  principal  requirement  of  any  (quasi-)  one-dimensional  solution  is 
to  quantize  the  accuracy  and  convergence  performance  of  the  finite  element 
algorithm.  Mixed  subsonic-supersonic  flows  are  particularly  useful,  since 
shock  sharpness  and  associated  undershoot/overshoot  are  readily  assessable. 
The  Riemann  shock  tube  problem,  cf.  [30,  p.  1007],  is  well  suited  to  this 
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analysis,  since  the  resultant  flow  structure  is  richly  endowed  with  discon¬ 
tinuities  and  sharp  variable  gradients,  interspersed  with  planar  plateau 
regions.  The  interesting  case,  with  unique  stagnation  sound  speeds  in  the 
two  chambers  separated  by  the  diaphragm,  has  been  exhaustively  examined  for 
algorithm  performance  in  the  finite  difference  field,  cf.  Van  Lear  [31]  and 
Zalesak  [32].  In  particular,  Sod  [33]  compares  the  results  of  approximately 
a  dozen  finite  difference  algorithms,  based  upon  Lagrangian  and  mixed  Lagran- 
gian-Eulerian  frameworks,  for  one  shock  tube  specification. 

The  developed  finite  element  algorithm  is  strictly  Eulerian,  and  there¬ 
fore  not  at  all  fine-tuned  for  the  non-linear  hyperbolic  inviscid  Riemann 
problem.  Therefore,  this  problem  serves  the  required  purpose,  as  well  as 
affording  an  interesting  exact  solution  for  comparison  to  assist  in  refining 
the  dissipation  parameter  set  The  diaphragm  is  placed  mid-way  in  a  duct 
of  uniform  cross-section,  with  the  unit  length  discretized  into  12  <  M  <  400 
finite  elements  R*  of  uniform  measure  A.  The  k=l,  M=100  discretization  corre¬ 
sponds  to  the  case  of  Sod  [33],  with  the  initial  condition  specifications 
u(x)  =  0,  p  =  1  =  p  on  0  <  x  <  0.5,  p  =  0.1  and  p  =  0.125  on  0.5  <  x  <  1.0 
with  y  =  1.4.  Figure  4  graphs  the  k=l  solution  variable  set  (QI(nAt)},  at 
t  =  0.14154s,  as  obtained  using  the  order-of-accuracy  optimized  parameter  set 

v1  s.  v  £  v2 ,  v1  =  0, v=  (15)~2.  Each  symbol  corresponds  to  a  nodal  coordinate  of 

q,  ex  m 

q  ,  and  the  dashed  lines  denote  the  initial  conditions.  The  shock  is  centered 

a 

at  x  =  0.75,  the  contact  discontinuity  is  centered  at  x  =  0.62,  and  the  rare¬ 
faction  wave  lies  upstream  of  x  =  0.5.  For  comparison.  Figure  5  is  a  graph 
of  the  M=99  solution  for  pn  and  p  generated  by  the  Lagrangian-rezone  Eulerian 
"MUSCL"  finite  difference  algorithm  of  Van  Lear  [31],  compared  to  the  analytical 
solution  shown  as  solid  lines.  The  k=l  finite  element  solution  is  considerably 
"smoother,"  with  less  well  defined  gradients  and  plateaus,  and  with  excessive 
overshoot/undershoot  about  the  shock.  Figure  6  is  the  improved  k=l,  M=200 
finite  element  algorithm  solution,  as  obtained  with  the  "numerically  optimized" 
dissipation  parameter  set  v1  =  v{3/8, 0,1/4}  and  v2  =  v{3/4,2,l},  for 

J.  ot  ok 

v  =  (15)  2  and  1  <  a  <  3.  Each  of  the  characteristic  Riemann  solution  features 
is  quite  accurately  represented  including  accurately  planar  plateaus  and  a 
crisply  defined  shock  with  negligible  overshoot. 
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Figure  5.  Solution  for  Riemann  Shock  Tube  Generated  By  the  MUSCL  Code, 
Reported  by  Van  iLeer  [31],  Courant  No.  =  0.9,  Ax  ='0;01, 
t  =  0.14154s. 


A  definitive  accuracy  assessment  can  be  made  in  terms  of  the  velocity 
and  internal  energy  distributions  computed  from  the  solution  set.  Figure  7 
compares  the  M=99  MUSCL  code  and  analytical  solution  (solid  line.  Figure  7a), 
with  the  va  optimized  k=l,  M=200  and  M=100  finite  element  solutions.  Figure 
7b)-7c)  respectively.  In  addition.  Figure  7d)  shows  the  M=100  data  obtained 
using  the  Crank-Nicolson  finite  difference  algorithm  equivalent  of  the  k=l 
finite  element  algorithm.  This  familiar  finite  difference  construction  as 
obtained  by  diagonalizing  the  [A200]  matrices  in  equation  107-109  and  111- 
113,  which  also  results  in  identical  elimination  of  the  v1  terms,  cf.  [23]. 
The  illustrated  solution  was  generated  using  v*  =  increasing  the 

dissipation  levels  would  moderate  the  illustrated  excessive  shock  overshoot, 
at  the  expense  of  farther  diffusing  the  shock  over  more  than  the  present  five 
domains.  In  comparison,  the  M=100,  k=l  finite  element  algorithm  results 
interpolate  the  shock  over  only  two  element  domains  with  negligible  over¬ 
shoot,  for  the  identical  CPU  and  main  memory  machine  requirements.  The 
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/ 


Figure  6. 


M=200,  k=l  Finite  Element  Solution, 
v£  =  v{3/8, 0,1/4},  v*  =  v{3/4,2 ,1} , 
Denotes  Initial  Conditions. 


Riemann  Shock  Tube 
t  =  0.14154s  (— ) 
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b)  M=200,  k=l  Finite  Element  Solution,  v 
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opt' 


Figure  7.  Finite  Element  and  Finite  Difference  Algorithm  Solution  Com¬ 
parisons,  Riemann  Shock  Tube,  t  =  0.14154s. 


48 


d)  Crank-Nicolson  Finite  Difference  Solution,  v1  =  0.,  v2  =  (15)”*. 

a  a  '  ' 


Figure  7.  Finite  Element  and  Finite  Difference  Algorithm  Solution  Com¬ 
parisons,  Riemann  Shock  Tube,  t  =  0.14154s,  Concluded. 
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M=100,  k=l  finite  element  data  compare  quite  favorably  to  the  MUSCL  code 
prediction;  the  velocity  resolution  of  the  shock  and  the  planarity  in  the 
high  temperature  plateau  are  nominally  identical.  The  MUSCL  code  data 
interpolates  the  contact  discontinuity  over  five  domains,  while  the  k=l 
data  smeared  it  over  nine  finite  elements.  Considering  that  the  Eulerian 
finite  element  algorithm  is  not  at  all  "hard-wired"  for  this  problem, 
the  evidenced  accuracy.is  excellent. 

The  performance  of  the  k=2  finite  element  algorithm  for  the  Riemann 
problem  is  nominally  comparable.  Following  numerical  experimentation,  the 
"optimized"  dissipation  parameter  set  was  determined  as  v1  =  0.  and  v2  = 

CX 

v{ 1/4 ,3/4,1/21.  Figure  8  graphs  the  resulting  solution  (QI(nAt)},  at 
t  =  0.14154s,  as  obtained  on  the  M=50  discretization  which  contains  as  many 
nodal  coordinates  as  the  k=l,  M=100  discretization.  Figure  9a  graphs 
the  resulting  solution  parameters  velocity  and  internal  energy,  and  the 
solid  lines  are  traces  of  the  exact  solution.  In  comparison  to  the  data 
of  Figure  7a)-c),  accuracy  is  excellent,  with  an  improvement  in  definition 
of  the  contact  discontinuity  evidenced.  In  the  case  of  the  k=2  algorithm, 
and  since  v^=0  is  "optimal,"  the  finite  difference  equivalent  diagonalized 
algorithm  results  are  not  severely  degraded.  Figure  9b). 

The  finite  element  algorithm  Jacobian  construction,  equations 
111-114,  exhibits  the  favorable  quadratic  convergence  rate  associated 
with  Newton  iteration.  For  example,  Table  2  presents  extremum  elements  of 
{ 6QI }  for  a  typical  integration  step  involving  three  iterations,  and  the 
nodal  location  of  these  extrema.  Convergence  is  at  least  quadratic,  including 
the  algebraic  pressure  equation.  The  sharply  defined  solutions,  as  obtained 
using  minimal  levels  for  v2,  typically  require  250  iterations  to  reach 
t  =  0.14154s,  at  a  Courant  number  (C  =  (u+a| At/Ax)  of  approximately  0.33. 

Since  the  algorithm  is  implicit,  the  sole  constraint  on  Courant  number 
is  accuracy.  For  the  discussed  results,  generated  for  C  =  0.33,  a  Richardson 
extrapolation  step  at  t  =  0.14154s  confirmed  that  the  significant  digit  in 

L 

||qn||  was  unaffected  by  the  temporal  truncation  error  associated  with  C,  as 
required  to  estimate  the  semi-discrete  (finite  element)  approximation  error. 
Figure  10  graphs  the  semi -di screte  convergence  characteristics,  measured  in 
both  H1'  and  E,  of  the  k=l  and  k=2  algorithm  solutions  for  the  Riemann  prob¬ 
lem,  over  a  discretization  range. 
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Figure  8  .  M=50,  k=2  Finite 

Shock  Tube,  t  =  0.14154 


X  COORDINATE 


Element  Algorithm  Solution,  Riemann 
,  v*  =  0.,  v*  =  v{l/4,  3/4,  1/2}. 
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b)..  Diagonalized  Quadratic,  v1  =  0.,  v2  =  v{l/4,  3/4,  1/2} 
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Figure  9.  Finite  Element  and  Diagonalized  Finite  Element  Algorithm  Com 
pari  sons,  Riemann  Shock  Tube,  M=50,  k=2,  t  =  0.14154s. 


TABLE  2 


NEWTON  ITERATION  CONVERGENCE  IN  {6QI} 
RIEMANN  SHOCK  TUBE 
M=100 ,  k=l,  e  =  O.OOl 


Iteration 

{6R}max 

Node 

(6M}max 

Node 

1 

-0.87325E-02 

70 

-0. 20780E-01 

70 

2 

0. 13072E-02 

69 

0 . 31246E-02 

69 

3 

0 . 2 1915E-03 

69 

0. 69326E-03 

69 

Iteration 

(6G}max 

Node 

{6P}max 

Node 

1 

-0 . 23351E-01 

70 

-0. 75154E-01 

70 

2 

0.34658E-02 

70 

-0 . 51826E-02 

71 

3 

0.46731E-03 

70 

-0 . 50114E-03 

68 

The  dashed  1 ines,  interpolating  the  k=2  data,  lie  uniformly  above  the  k=l 
data  (solid  lines),  for  each  dependent  variable,  and  each  has  a  nominally 
identical  slope.  A  modest  slope  distribution  exists  among  the  various 
dependent  variable  data.  Assuming  this  distinction  insignificant,  a  nominal 
convergence  rate  of  2/3  is  indicated,  in  qualitative  agreement  with  a  dif- 

X  L. 

ference  scheme  of  p  '  order  accuracy,  ie., 

II*  <  ci  ♦  ...  (115) 

where  p  =  2  and  Ci  is  a  constant  independent  of  A,  the  measure  of  the  uniform 
mesh.  The  three  solid  symbols  on  A0  =  10"2  are  the  data  in  .E  generated 
by  the  M=100  finite  difference  solution  shown  in  Figure  7d).  The  consistent 
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©v  r SI  ope  =  2/3 


trend  in  error  extremi zation  by  a  finite  element  algorithm  thus  appears  con¬ 
firmed  for  this  problem , (cl ass ) . 

Solution  efficiency  is  principally  determined  by  discretization  (M)  and 
integration  timestep  (C),  in  concert  with  the  requested  convergence  level 
for  the  Newton  iteration  algorithm.  Table  3  summarizes  the  effect  of 
decreasing  solution  cost  by  increasing  C  for  the  M=100,  k=l  Riemann  solution. 
Solution  cost  can  be  nominally  halved  by  increasing  the  Courant  number  to  a 
certain  point.  Thereafter,  diminishing  returns  are  encountered  as  additional 
Newton  iterations  are  required  for  convergence.  In  all  data  for  C  >  0.33, 
the  significant  digit  in  the  semi -discrete  approximation  norm  is  being  directly 
affected  by  temporal  truncation  error,  and  solution  accuracy  is  degraded 
monotonically  with  increasing  C.  Of  particular  note,  the  k=l  finite  element 
algorithm  solution  accuracy  at  C  =  1.41  is  uniformly  superior  to  the  C=0.33 
finite  difference  solution,  see  Figure  10,  and  was  obtained  at  40%  of  the  com¬ 
puter  cost. 


TABLE  3 

ACCURACY  AND  EFFICIENCY  SUMMARY 
RIEMANN  SHOCK  TUBE 
M=100 ,  k=l,  e=0.001 


Courant 

Number 

No.  Time 
Steps 

No.  of 
Algorithm 
Passes 

GO  Step 
CPU 

formalized) 

L. 

Energy  Norm  ||q"||E. 

P 

m 

g 

0.33 

80 

229 

1.00 

1  :  1 

2.43 

17.9 

0.44 

60 

180 

0.79 

2.39 

17.4 

0.67 

40 

121 

0.54 

tlllll 

2.29 

16.8 

0.90 

30 

92 

0.42 

1.41 

2.15 

16.0 

1.41 

20 

82 

0.38 

1.34 

1.90 

14.4 
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3.  Mixed  Flow  In  A  Variable  Cross-Section  Duct 


Since  the  finite  element  algorithm  is  implicit,  all  solution  zones 
are  in  continuous  communication,  even  those  regions  where  the  flow  is 
locally  supersonic.  An  interesting  test  occurs  for  a  mixed  flow  situation, 
wherein  subsonic  boundary  data  modification  require  an  adjustment  to  the  flow 
in  a  supersonic  region.  As  an  example,  off-design  operation  of  a  de  Laval 
nozzle  yields  a  strong  shock  in  the  divergent  section,  such  that  the  shocked 
flow  can  diffuse  to  the  exit  chamber  pressure.  A  modest  increase  in  the 
chamber  pressure  requires  the  shock  to  relocate  upstream  of  its  original 
location,  wherein  the  flow  was  supersonic,  with  the  shock  always  separating 
the  mixed  flow  regions.  Anytime  during  this  period  of  adjustment,  the  super¬ 
sonic  flow  upstream  of  the  shock  must  not  respond  to  the  new  chamber  pressure 
level  ,  the  shock  must  progressively  weaken  as  it  moves  upstream  through 
regions  of  smaller  cross-section,  and  the  downstream  flow  must  diffuse 
smoothly  to  the  new  exit  conditions. 

Figure  11  summarizes  the  M=74,  k=l  algorithm  steady-state  solution, 
for  an  off-design  nozzle  flow  with  Mi  =  1.35  upstream  of  the  shock,  using 
vo^  and  i  v  The  solid  line  is  the  exact  solution;  the  1  solution  is 
in  excellent  agreement  on  shock  strength  and  definition.  It  is  interesting  how 
the  momentum  solution  attempts  to  interpolate  the  discontinuity  of  the  shock, 
which  of  course  will  be  of  finite  width.  The  subsonic  inlet  boundary  conditions 
are  specified  p,  m,  and  g  with  p  computed  from  equation  7.  The  subsonic  out¬ 
let  condition  is  p  specified,  and  a  vanishing  normal  derivative  for  p,  m,  and  g. 

Using  the  solution  of  Figure  11  as  Initial  conditions,  the  specified 
exit  pressure  was  raised  \b%  and  held  fixed  at  the  new  level.  The  inlet  condi¬ 
tions  remained  as  originally  specified.  Figure  12  graphs  the  k=l  algorithm  solu¬ 
tion  at  the  new  steady-state,  where  the  solid  lines  are  the  initial  conditions. 

The  flow  has  adjusted  fully  to  the  new  exit  pressure,  with  a  new  shock  Mach 
number  of  Mi  =  1.15.  The  flow  upstream  of  the  shock  is  unaltered  from  the  initial 
conditions,  while  a  smooth  subsonic  diffusion  to  exit  conditions  is  predicted 
downstream.  The  shock  is  interpolated  across  four  elements  with  no  overshoot. 
Figure  13  summarizes  the  comparison  M=37,  k=2  algorithm  solution,  as  obtained 
using  =  £v  £.  The  solution  accuracies  are  indistinguishable. 
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X  COORDINATE 


Figure  12.  M=74,  k=l  Algorithm  Solution  for  Mixed  Flow  in  a  Variable 

Cross-Section  Duct,  \>a  =  ivQpt >  (  )  Initial  Condition. 
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X  COCSMNJOE  X  COORDINATE 

Figure  13.  M=37,  k=2  Algorithm  Solution  for  Mixed  Flow  in  a  Variable 

Cross-Section  Duct,  ( - )  Initial  Condition. 
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4.  Two-Dimensional  Flow,  Metric  Data 


As  noted  in  the  Introduction,  a  wide  variety  of  methodologies  exist 

A 

for  generation  of  the  inverse  coordinate  transformation  n-  =  n-(x-),  see 

1  J 

equation  28.  For  example.  Figure  14  illustrates  body-fitted  coordinate 
systems  for  various  two-dimensional  aerodynamic  configurations,  as  generated 
using  Poisson  equation  solutions  [8].  The  principal  requirement  is  to  make 
contours  corresponding  to  aerodynamic  surfaces,  and/or  symmetry  or  periodocity 
boundaries,  be  coordinate  curves  of  the  system.  The  second  requirement, 
from  the  standpoint  of  simulation  accuracy,  is  to  avoid  highly  distorted 
computational  cells  in  critical  flow  regions,  eg.,  stagnation  points. 

With  respect  to  the  generalized  coordinates  finite  element  algorithm, 

each  of  the  grids  illustrated  in  Figure  14  is  the  definition  of  the  array 

of  nodal  coordinates  {XI}  =  Z {XI }  ,  l<e<M,  1<I<2,  required  to 

e  e  -  -  - 

define  the  coordinate  transformation,  see  equations  29-32.  Restricting 
attention  to  the  bilinear  basis  { N i ( n ) > ,  k=l  in  equation  29,  the  typical 
finite  element  domain  R|  contains  only  vertex  nodes.  For  the  counterclock¬ 
wise  numbering  convention,  Figure  1,  denote  the  element  node  coordinates  as 
{XII  =  {XI,  YI,  1  <  I  <  4}.  The  array  {ETAKJL,  1  <  (K,J)  <  2,  is  con- 

c  —  —  c  —  — 

structed  via  equation  34,  and  the  calculus  operation  is  elementary.  Define 
the  discrete  approximation  on  R|  using  {Ni(n)>  as 


fs%l 


E  ddfj  {Ni>T{ETAKJ}e,  1  <  (K,J)  <  2 


(116) 


Noting  that  (det  J)”1  cancels  in  the  algebra,  the  elements  of  {ETAKJ}g  are 


{ETAKJ}g  =  \ 


{Y4-Y1 ,  Y3-Y2 ,  Y3-Y2,  Y4-Y1}  , 
{X1-X4,  X2-X3,  X2-X3,  Xl-X4}„, 
{Y1-Y2 ,  Y1-Y2  ,  Y4-Y3,  Y4-Y3F, 
{X2-X1 ,  X2-X1 ,  X3-X4,  X3-X4}°, 

. 


(1,1) 

(1,2) 

(2,1) 

(2,2) 


' 

?  (117) 


where  the  numbers  in  parenthesis  indicate  (K,J). 
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b)  Highly  Cambered  Airfoil  and  Trailing  Edge  Close-Up,  from  Sorenson 
and  Steger  [8,  p.  456] 

Figure  14.  Examples  of  Computational  Grid  Transformations  For  Two- 
Dimensional  Aerodynamic  Flow  Predictions. 
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c)  C-Type  Grid  For  Airfoil  in  a  Wind  Tunnel,  from  Sorenson  and  Steger 
[8,  p.  458] 


d)  Leading  Edge  Detail  e)  C-Type  Grid  For  Nacelle  With  Center- 

Sorenson  and  Steger  [8,  p.  457]  body,  from  Kowalski  [8,  p.  349]. 


Figure  14.  Examples  of  Computational  Grid  Transformation  For  Two-Dimen¬ 
sional  Aerodynamic  Flow  Prediction,  Concluded. 
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The  corresponding  definition  for  det[J]g  is 


det[J]e  =  {N1}T{DET}e 


(118) 


and 


(X2-X1) (Y4-Y1) 
J (X2-X1) (Y3-Y2) 
(X3-X4) (Y3-Y2) 
(X3-X4) (Y4-Y1) 


(X4-X1) (Y2-Y1) 
(X3-X2) (Y2-Y1) l 
(X3-X2) (Y3-Y4) 
(X4-X1) (Y3-Y4) 


(119) 


The 

on  R2  is, 
e 


r.h 


eontrava.riant  convection  velocity  approximation  uj^ 


equation  39, 


uk6  5  {N!}T{UBARK}e  (120) 

and  the  algebra  operations  yield. 

{UBARK)e  =  {ETAKJ>e  TNi>T{UJ>e  (121) 

Equation  121  is  evaluated  on  a  nodal  basis,  whereupon  {Nil  reduces  to  the 
Kronecker  delta.  Summation  is  implied  over  J,  and  at  the  nodes  UJ  =  MJ/R, 
ie. ,  u^  =  rn./p. 

The  tensor  product  Jacobian  formulation,  given  for  the  quasi-one-dimen- 
sional  flow,  is  modified  only  to  account  for  alignment  with  either  the  ni  or 
n2  coordinate  axis.  The  cardinal  basis  remains  as  given  in  equation  104,  and 
DETg  is  equal  to  one-half  the  element  measure,  see  equation  106.  Hence,  for 
equations  66-71,  and  using  Pythagorus1  rule 


(122) 

where  R  and  L  denote  "right"  and  "left"  respectively,  for  realignment.  Using 
the  same  considerations,  the  elements  of  {ETAKJ)e  are  equivalent  to  the  direction 
cosines  of  the  local  affine  transformation  between  x.  and  q..  Denoting  Gas 

*  J 


{DETKL  ^  i/(XKR  -  XKL)2  +  (YKR  -  YKL): 
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the  angle  between  the  xx  and  nx  axes  on  R2,  for  example,  then 


{ETAKJlg 


{  cos0,  sin©},  (1,1) 
{-sin©,  cos©},  (1,2) 


(123) 


Using  equation  123,  the  elemental  definition  for  convection  velocity  in 
[Jn]  remains 


{UBARK>e  =  {ETAKJ}e{Ni)T{UJ}e  (124) 

5.  Phase  Accuracy,  The  Rotating  Cone 

The  rotating  cone  test  case  is  a  multi -dimensional  solution  of  the 
linear  hyperbolic  continuity  equation  1,  A  given  initial  density  distribu¬ 
tion  p°(x,0)  is  subject  to  an  imposed  constant  convection  velocity  field 
U0  =  raj©,  where  w  is  the  angular  velocity.  The  analytical  solution  is  exact 
preservation  of  p°  for  all  time,  and  the  test  is  extremely  demanding  with 
respect  to  algorithm  phase  accuracy  and  embedded  artificial  diffusion. 

The  initial  distribution  will  be  rapidly  destroyed  by  either  error  mechanism. 

The  finite  element  tensor  product  algorithm,  with  zero  artificial  dissi 
pation,  e  (J  in  equation  22,  can  produce  solutions  to  the  rotating  cone 
test  problem  exhibiting  excellent  accuracy.  The  solution  domain  R2  is  the 
unit  square,  discretized  uniformly  into  an  M  =  32  x  32  mesh,  for  the  k=l 
algorithm,  and  an  M=16xl6  mesh  for  the  k=2  algorithm.  Both  meshes  thus 
possess  33x33  =  1089  nodal  coordinates.  The  initial  density  distribution 
is  a  "cosine  hill,"  obtained  by  rotating  a  cosine  distribution  about  the  peak 
level.  Figure  15a  illustrates  the  resultant  distribution  for  p°(x,y),  as 
obtained  for  interpolating  the  cosine  half-wave  onto  seven  nodal  coordinates 
(with  levels  15,  50,  85,  100,  85,  50,  15).  A  non-reflecting,  vanishing  nor¬ 
mal  derivative  boundary  condition  is  employed,  a^  =  0  =  a*  in  equation  18. 

Figure  15  summarizes  the  various  numerical  solutions,  as  a  function 
of  Courant  number  C,  and  semi-discrete  approximation  function  degree  k,  fol¬ 
lowing  2t r  rotation  from  the  initial  condition  specification.  For  comparison, 
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Figure  15.  Rotational  Convection  of  a  Cosine  Hill  Distribution. 


65 


Figure  15a)  is  also  the  exact  analytical  solution,  and  shows  the  sense  of 
rotation  induced  by  U  .  Figure  15b)  is  the  solution  generated  by  the  non- 
dissipative  k=2  algorithm,  using  a  constant  integration  time  step  At 
yielding  C  =  0.25  at  the  center  of  p(x,t).  (Note  that  C  varies  linearly, 
from  zero  at  the  center  of  R2,  as  a  consequence  of  U0.)  The  accuracy  of  this 
solution  is  excellent;  the  final  distribution  exhibits  almost  exactly  the 
symmetries  of  p° ,  the  peak  level  has  actually  increased  slightly  to  105, 
and  the  extremum  magnitude  of  the  trailing  phase  dispersion  error  wake  is 
2%  and  quantized  by  the  height  of  the  plot  base  of  R2.  The  modest  back¬ 
ground  waviness,  on  the  plane  exterior  to  p°,  is  the  residual  dispersion 
induced  by  the  passage  of  the  solution. 

Increasing  the  integration  time-step  At(C),  and/or  decreasing  the 
semi -discrete  approximation  degree  k,  degrades  the  accuracy  of  the  solution. 

For  example.  Figure  15c)  graphs  the  final  solution  p(x,t^),  as  obtained  by 
the  k=2  algorithm  at  C  =  0.5,  and/or  the  k=l  algorithm  at  C  =  0.25.  In  both 
cases,  the  additional  phase  distortion  has  retarded  the  solution  group  velocity, 
such  that  the  peak  has  not  reached  the  correct  nodal  coordinate  by  tf,  hence 
the  flattened  appearance.  However,  marching  forward  an  additional  one  or  two 
At  recovers  the  solution  extremum  of  100,  for  k=2,  and  98  for  k=l,  at  a  nodal 
coordinate.  The  maximum  amplitude  of  the  trailing  dispersion  wake  is  91, 
note  the  higher  plot  base,  and  the  background  ripple  structure  throughout  R2 
is  modestly  more  pronounced.  Figure  15d)  graphs  the  k=l  algorithm  solution 
obtained  for  C  =  0.5.  The  indicated  solution  peak  is  93,  and  is  located  one 
nodal  coordinate  shy  of  where  it  should  be.  The  extremum  dispersion  wake  error 
is  17%,  hence  the  rather  pronounced  plot  base.  However,  the  solution  accuracy 
is  still  quite  acceptable. 

For  comparison,  Figure  15e)  graphs  the  solution  obtained  for  C  =  0.25 
by  the  Crank-Nicolson  finite  difference  equivalent  of  the  non-dissipative  k=l 
algorithm.  This  algorithm  is  obtained  simply  by  diagonalizing  the  matrix 
[B200] ,  equation  74.  This  ad  hoc  operation  destroys  the  phase  accuracy,  recall 
Figure  3,  and  hence  yields  the  extremely  poor  solution  for  the  rotating  cone. 
The  peak  has  decayed  to  47,  the  extremum  wake  error  magnitude  is  25%  of  p“ 
and  the  phase  error  is  well  dispersed  over  R2.  For  completeness.  Figure  15f) 
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graphs  the  solution  for  the  Crank-Nicolson  equivalent  of  the  k=2  algorithm 
at  C  =  0.5.  The  peak  is  85,  the  wake  extremum  is  21%,  and  the  dispersion 
error  is  well  dispersed  onto  Rz. 

These  computational  results  confirm  the  phase  accuracy  robustness 
of  the  developed  finite  element  algorithm  construction.  This  attribute  is 
deemed  particularly  important  when  addressing  high  Reynolds  number  flow 
prediction.  The  results  of  several  additional  tests  are  reported  [15], 
wherein  various  combinations  of  v*  and  vz  were  employed  for  k=l  solutions  of 
the  rotating  cone  test.  The  magnitude  of  the  extremum  wake  error  could  be 
reduced,  using  the  dissipation  mechanism,  but  always  at  the  expense  of  dif¬ 
fusing  the  solution  peak  as  well.  From  the  standpoint  of  accuracy,  these 
data  strongly  confirm  that  poor  algorithm  phase  accuracy  cannot  be  corrected 
through  use  of  artificial  dissipation. 

6.  Two-Dimensional  Riemann  Shock  Tube. 

Several  additional  aspects  regarding  algorithm  accuracy  and  performance 
aspects  can  be  assessed  using  as  a  test  case  the  two-dimensional  equivalent  of 
the  Riemann  shock  tube.  Using  the  convergence  data  discussed  in  Section  V.2, 
a  reasonable  (affordable)  discretization  of  the  Riemann  problem  was  defined 
as  M=32xN,  where  4  <  N  <  20  dependent  upon  whether  the  test  problem  was 
inviscid  or  viscous.  For  reference,  Figure  16  graphs  the  one-dimensional, 
k=l  algorithm  solution  parameters  of  velocity  and  internal  energy  for  the 
Riemann  problem.  The  solid  lines  are  traces  of  the  M=200,  k=l  solution  from 
Figure  7b),  and  each  open  symbol  is  a  nodal  solution  value  obtained  using 
=  v^.  Interestingly,  on  the  progression  to  quite  coarse  grids,  the  impor¬ 
tance  of  vS  0  on  accuracy  becomes  diminished.  The  solid  symbols  in  Figure 
16  are  the  appropriate  nodal  values,  obtained  using  v*  =  0,  vZ  =  v{3/4,2,l) 
in  the  k=l  algorithm  solution,  at  locations  where  the  two  solutions  differed. 
As  can  be  seen,  the  differences  are  truly  negligible;  hence,  the  two-dimen¬ 
sional  Riemann  solutions  were  all  executed  using  v1  =  0  and  v2  =  v2  .  . 

a  a  opt 

Several  tests  were  conducted  to  evaluate  application  of  vanishing 
gradient  boundary  conditions,  on  the  transverse  walls  of  the  two-dimensional, 
inviscid  flow  Riemann  shock  tube.  Figure  17  graphs  the  M=32x6,  k=l  solution 
at  t  =  0.14154s,  for  the  shock  tube  aligned  parallel  to  the  Xi-coordinate 
direction.  As  a  consequence,  the  entire  momentum  solution  is  carried  by 
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a)  Density  p  b)  Momentum  mi  c)  Energy  g  d)  Pressure  p 

Figure  17.  M=32x6,  k=l  Finite  Element  Solution,  Two-Dimensional  Riemann 

Shock  Tube,  t  =  0.14154s,  v1  =  0,  v?  =  v{3/4,  2,  1}. 


mi i  indeed,  m2  was  computed  to  equal  zero  to  three  significant  digits,  for 
the  Newton  convergence  requirement  set  at  e  =  0.01.  The  M=32x6  solution  ■; 
agrees  almost  exactly  with  the  companion  M=32xl,  one-dimensional  k=l  algo¬ 
rithm  solution,  see  Figure  18,  even  to  prediction  of  the  small  amplitude 
"2Ax"  waves  in  the  zero  velocity  region  upstream  of  the  rarefaction  wave. 

As  shown  in  Figure  17,  the  tensor  product  algorithm  can  readily  enforce 
a  vanishing  derivative  transverse  wall  boundary  condition  without  use  of  ex¬ 
traneous  phantom  cells  exterior  to  R2.  A  more  demanding  case  corresponds 
to  insertion  of  the  shock  tube  diaphragm  at  an  angle,  which  yields  the 
requirement  for  the  algorithm  to  predict  a  shock  oblique  to  the  mesh.  The 
corresponding  M=32x6,  k=l  algorithm  solution  at  t  =  0.14154s  is  graphed  in 
Figure  19.  The  algorithm  again  exhibits  excellent  solution  fidelity,  and 
experienced  no  difficulty  in  enforcing  the  gradient  boundary  constraint  for 
steep  solution  gradients  oblique  to  the  walls.  There  occurs  a  slightly 
enhanced  "2Ax"  trashiness  in  the  solution  upper  right  corner,  and  a  modest 
no-zero  distribution  for  m2  is  also  computed. 

Figure  20  graphs  the  M=32x6,  k=l  algorithm  solution  for  the  situation 
where  the  shock  tube  is  aligned  at  an  angle  in  the  laboratory  reference 
frame.  As  a  consequence,  since  m^  is  resolved  into  scalar  components  parallel 
to  x.j ,  both  momentum  equations  generate  non-zero  solutions.  As  can  be  seen, 
comparing  Figures  17b)  and  20b)-20c),  the  shock  does  not  appear  too  notice¬ 
able  in  the  m^  solutions.  However,  computation  of  the  component  of  u^  parallel 
to  ni,  which  is  the  coordinate  parallel  to  the  shock  tube  axis,  using  equa¬ 
tion  121  with  {MJ}  and  {R},  confirms  the  existence  of  the  shock.  Figure  20e). 
Further,  the  prediction  of  the  internal  energy  distribution.  Figure  20f),  is 
in  good  agreement  with  the  "correct"  solution,  see  Figure  16. 

As  the  final  test,  the  two-dimensional,  aligned  shock  tube  definition 
of  Figure  17  was  repeated  as  a  viscous  problem  at  Re  s  105,  for  an  M=32x20 
uniform  discretization.  The  corresponding  no-slip  boundary  conditions  on  the 
transverse  walls  are  mi  =  0  =  m2.  The  cold  wall  boundary  condition  e  =  e0, 
equation  6,  was  defined  for  the  energy  equation,  and  a  vanishing  normal  deri¬ 
vative  applied  for  density.  No  boundary  condition  specifications  are  appro¬ 
priate  for  p,  a.,  or  q.,  since  each  is  defined  as  an  algebraic  equation.  The 
integration  time  step  and  Newton  convergence  requirement  were  maintained 
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X  CCOTOBttTE 


X  COOSOIKATt 


Figure  18.  M=32xl ,  k=l  Finite  Element  Solution  For  Riemann  Shock  Tube, 
( - )  M=200  Solution,  (•)  vj  =  0,  vj  =  v  t,  t  =  0.14154s. 
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a)  Momentum  mx 


d)  Energy  g 

c)  Pressure  p 

b)  Momentum  m2 


Figure  19.  M=32x6,  k=l  Finite  Element  Algorithm  Solution,  Two-Dimensional 

Riemann  Shock  Tube,  Oblique  Diaphragm,  t  =  0.14154s. 

identical  to  the  inviscid  test  case  values.  Figure  21  graphs  the  resulting 
k=l  algorithm  solution  at  t  =  0.14154s.  In  the  centroidal  region,  away  from 
the  influence  of  the  no-slip  wall,  the  viscous  solution  agrees  almost  exactly 
with  the  M=32,  k=l  one-dimensional  inviscid  solution.  Figure  18.  The  influ¬ 
ence  of  the  viscous  flow  boundary  conditions  is  clearly  evident  in  all  depen¬ 
dent  variables.  The  growth  of  a  viscous  boundary  layer  behind  the  travelling 
shock  is  just  visible  in  mx.  The  solution  for  m2  is  highly  oscillatory,  but 
the  peak  value  is  only  about  10%  of  the  maximum  mi.  The  o12  shear  stress 
solution  is  sharply  peaked  at  the  wall,  with  a  steep  front  adjacent  to  the 
shock.  The  solution  is  non-zero  only  in  the  vicinity  of  the  wall,  and  ex¬ 
hibits  the  required  skew-symmetries. 
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a)  Density  p  b)  Energy 


c)  Momentum  mi  d)  Momentum  m2 

Figure  21.  M=32x20,  k=l  Finite  Element  Solution,  Two-Dimensional  Viscous 

Riemann  Shock  Tube,  Re  =  105,  t  =  0.14154s. 
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e)  Pressure  p  f)  Shear  Stress  c12 


Figure  21.  M=32x20,  k=l  Finite  Element  Solution,  Two-Dimensional  Viscous 

Riemann  Shock  Tube,  Re  =  105,  t  =  0.14154s,  Concluded. 

This  test  problem  specification  possessed  ten  dependent  variables 
per  nodal  coordinate,  plus  five  pieces  of  metric  data,  yielding  approxi¬ 
mately  10,000  nodal  degrees  of  freedom  to  handle.  The  Newton  convergence 
character  was  basically  identical  to  the  inviscid  flow  test  case  experience, 
following  a  few  extra  iterations  at  start-up.  The  solution  employed  27 
time  steps  to  reach  t^,  required  90  iterations,  and  utilized  70K  words  of 
central  memory. 
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7.  Three-Dimensional  Problems,  Metric  Data 

The  metric  data  for  three-dimensional  problem  definitions  is  generated 
in  the  manner  discussed  for  R2  in  Section  V.4.  Denote  the  array  of  nodal 
coordinate  triples  of  the  discretization  UR|  as  {XI}  =  E  {XI }g ,  1  <  e  <  M 
and  1  <  I  <  3.  For  the  k=l  basis,  recall  Figure  lb),  the  finite  element 
domain  R|  possesses  only  the  vertex  nodes  {XI }g  =  {XI ,  YI ,  ZI ,  1  <  I  <  8}. 

For  the  definition  given  in  equation  116,  the  elements  on  the  trace  of 
[3iv/8x-],  evaluated  at  the  centroid  (ri .  =  0,  l<i  <3)  of  R2  are 

K  J  1  -  -  g 

{ETAll}^  =  HY41Z51,  Y32Z62 ,  Y32Z73,  Y41Z84, 
e  Y85Z51 ,  Y76Z62 ,  Y76Z73,  Y85Z84}e 

{ETA22 }T  =  HX21Z51,  X21Z62,  X34Z73,  X34Z84, 
e  X65Z51 ,  X65Z62 ,  X78Z73,  X78Z84>e 

{ETA33l[  =  HX21Y41,  X21Y32,  X34Y32,  X34Y41 

6  X65Y85,  X65Y76,  X78Y76,  X78Y85}e  (125) 


For  a  rectangular  domain,  {ETAKJ>e  ={0}  for  K  f  J,  and  the  notation  in 
equation  125  is  defined  as,  for  example 

Y41Z51  =  (Y4  -  Y1)(Z5  -  Zl)  (126) 


In  the  same  notation,  the  elements  of  {DET}g,  evaluated  at  the 
origin  of  the  domain  R|,  are 

{DET}1  =  i{X21Y41Z51,  X21Y32Z62,  X34Y32Z73,  X34Y41Z84, 

9  X65Y85Z51 ,  X65Y76Z62,  iX78Y76Z73,  X78Y85Z84>e  (127) 

In  the  instance  of  R?  being  a  rectangular  parallelepiped,  then  {DET1  = 
Ae{l}/64,  where  Ag  is  the  volume.  For  equations  125-127,  nodes  1-4  lie  in 
the  lower  plane,  see  Figure  lb,  with  node  1  in  the  lower  left  corner.  Nodes 
5-8  are  also  ordered  counterclockwise,  in  the  upper  plane,  with  node  5  above 
node  1. 
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For  the  tensor  product  Jacobians,  equations  66-71, 


{DETK)e  =  MXKR  -  XKL)2  +  (VKR  -  YKR)2  +  (ZKR  -  ZKL)2  jjj  (128) 

where  R  and  L  denote  "right"  and  "left"  respectively,  for  the  specific 
alignment.  As  for  R2,  the  components  of  (ETAKJ)e  for  equations  66-71  are 
equivalents  the  di  rection  cosines  of  the  local  affine  transformation  between  x- 
and  q.,  1  <  ( l , j )  <  3.  The  definition  for  (UBARK)  remains  unchanged  from 

J  "  ™  ® 

equation  124. 

8.  Three-Dimensional  Riemann  Shock  Tube 

The  developed  tensor  product  algorithm  has  been  evaluated  for  the 
three-dimensional  equivalent  of  the  inviscid  Riemann  shock-tube  problem, 
using  the  k=l  basis  on  an  M  =  32x4x4  uniform  discretization  of  R3.  The 
principal  requirement  is  to  quantize  performance  of  the  tensor  product 
Jacobian  and  the  various  vanishing  normal  gradient  boundary  conditions.  The 
shock, tube  is  aligned  with  the  xx  axis,  and  the  dissipation  parameter 
set  is  retained  as  defined  for  the  two-dimensional  simulation,  v*  =  0, 
v*  =  v{3/4,2 ,2 ,2,1}. 

Each  x3  plane  of  the  solution  exhibits  identical  solution  distributions, 
including  at  the  nodal  coordinates  defined  on  the  boundary  3R  of  R3.  Further, 
each  distribution  is  identical  with  the  comparison  two-dimensional  solution, 
see  Figure  17.  For  example,  Figure  22  shows  the  actual  print- out  for 
mi(xi,x2),  on  tihe  planes  x3  =  constant,  at  t  =  0.14154s.  The  total  number 
of  degrees  of  freedom  for  this  case  was  approximately  3200,  excluding  the 
metric  data.  The  solution  employed  27  steps  to  reach  t^,  required  56 
iterations  and  utilized  64K  words  of  central  memory. 
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Figure  22.  Program  Print-Out  For  Three-Dimensional  Riemann  Shock  Tube,  Principal 
Momentum  Component  mls  t  =  0.14154s. 
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SECTION  VI 


CONCLUSIONS 

An  implicit,  generalized-coordinates  finite  element  numerical  solution 
algorithm  has  been  derived  and  evaluated  for  the  three-dimensional  Navier- 
Stokes  equations.  The  theoretical  basis  utilizes  a  Galerkin-Weighted  Residuals 
statement,  rendering  the  semi-discrete  error  orthogonal  to  the  finite  element 
approximation  subspace,  augmented  with  a  penalty  constraint  forcing  orthogon¬ 
ality  of  the  gradient  of  this  error  as  well.  As  a  consequence,  the  algorithm 
possesses  highly  phase  selective  dissipation  mechanisms  permitting  accurate 
resolution  of  solutions  exhibiting  a  high  degree  of  non-smoothness.  A  Fourier 
stability  analysis  yields  an  estimate  of  the  dissipation  parameter  set, 
which  has  been  refined  to  enhance  accuracy  for  a  shocked  flow  prediction. 

Multiple  factors  affecting  eolation  accuracy,- convergence  and  effi¬ 
ciency  have  been  examined.  The  generalized  coordinates  framework  directly 
facilitates  matching  of  arbitrary  surface  descriptions  of  the  solution 
domain  for  complete  geometric  versatility.  This  formulation  as  well  permits 
establishment  of  the  tensor  matrix  product  approximation  to  the  Newton  itera¬ 
tion  algorithm  Jacobian,  reducing  by  orders  of  magnitude  the  memory  and  CPU 
requirements  for  a  multi-dimensional  problem  definition.  The  algebraic  and 
constitutive  equations,  defining  pressure,  stress  tensor  and  heat  flux  vector, 
are  handled  in  an  identical  manner,  yielding  an  overall  consistency  to  the 
algorithm. 

The  results  of  a  variety  of  numerical  experiments  have  verified  the 
basic  attributes  of  the  developed  algorithm.  However,  with  the  limited 
time  frame  and  funding  for  this  project,  the  scope  of  numerical  results  that 
could  be  assessed  has  barely  been  touched.  The  theoretical  construction  of 
the  algorithm  does  appear  to  contain,  as  a  special  case,  the  popular  approxi¬ 
mate-factorization  finite  difference  algorithms  in  use  today.  However,  the 
troublesome  viscosity  mixed  derivatives  are  absent  from  the  finite  element 
formulation,  yet  the  impact  of  this  dependent  variable  choice  has  not  been 
evaluated.  Several  other  comparisons  could  and  should  be  made,  to  more  firmly 
quantify  the  desired  algorithm  comparisons. 
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APPENDIX  A 


FINITE  ELEMENT  ALGORITHM  HYPERMATRICES 
LINEAR  AND  QUADRATIC  BASIS  ON  ONE  DIMENSIONAL  SPACE 

1.  Linear  Finite  Element  Formulation 
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2.  Quadratic  Finite  Element  Formulation 
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APPENDIX  B 


FINITE  ELEMENT  ALGORITHM  HYPERMATRICES 
BI-LINEAR  TENSOR  PRODUCT  BASIS  ON  TWO-DIMENSIONAL  SPACE 
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APPENDIX  C 


FINITE  ELEMENT  ALGORITHM  MATRICES 
TRI -LINEAR  TENSOR  PRODUCT  BASIS  ON  THREE-DIEMNSIONAL  SPACE 
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